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Quantum theory of the third-order nonlinear electrodynamic effects in graphene 

S. A. Mikhailo\0 

Institute of Physics, University of Augsburg, D-86135 Augsburg, Germany 
(Dated: February 3, 2016) 

The linear energy dispersion of graphene electrons leads to a strongly nonlinear electromagnetic 
response of this material. We develop a general quantum theory of the third-order nonlinear local 
dynamic conductivity of graphene aag~f5{Lij\,LiJ2,oJs), which describes its nonlinear response to a 
uniform electromagnetic field. The derived analytical formulas describe a large number of different 
nonlinear phenomena such as the third harmonic generation, the four wave mixing, the saturable 
absorption, the second harmonic generation stimulated by a dc electric current, etc., which may be 
used in different terahertz and optoelectronic devices. 
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I. INTRODUCTION 


Electrodynamic properties of graphene have attracted considerable attention in recent yearsi^— . Microwave and 
optical response, plasma and cyclotron resonances, nonlinear electromagnetic properties have been investigated both 
theoretically and experimentally in many research groups in the world^^— . Special interest to these topics is generated 
by the opportunity to use many of these fundamental phenomena for microwave-, terahertz- and optoelectronic 
applications^i^d^. 

The most distinctive feature of graphene is the linear energy dispersion of its quasi-particles, electrons and holes, 
E±{p) = ±vf\p\- It was theoretically predictedi^ and then experimentally confirme d^°’^^ , that due to this property 
graphene should demonstrate strongly nonlinear electromagnetic behavior. Physically this can be understood as 
followsiS. Assume that a particle with the linear spectrum is placed in the uniform external electric field E{t) = 
Eo cos Lut. According to Newton equations of motion its momentum will oscillate as p{t) oc sinwt, while the velocity 
v{t) = dE/dp oc p{t)/\p{t)\, as well as the current j{t), which are not proportional to the momentum, will be strongly 
nonlinear functions of p{t) and will therefore contain higher frequency harmonics, 


j{t) oc v{t) oc sgn(sinajt) oc sinwt -I- - sin3wt - 

o 


( 1 ) 


In other words, since graphene electrons have no mass and can move in any direction only with the velocity vp sa 10® 
cm/s, at the return points they have to instantaneously change their velocity from -1-10® cm/s to —10® cm/s which is 
accompanied by a strong acceleration and leads to emission of electromagnetic radiation. The nonlinear electrodynamic 
phenomena in graphene have been further studied in Refsi^^— (theory) and in Refsi^^— (experiment). 

In the first theoretical publications the nonlinear (third-order) electromagnetic response of graphene was studied in 
the quasiclassical approximation, using the Boltzmann kinetic approach. Such a theory takes into account only the 
mtra-band oscillations of the graphene electrons, which is valid at low (microwave/terahertz) frequencies hio < 2Ef, 
where w is the radiation frequency, Ep = \p\ is the Fermi energy and p, is the chemical potential. At higher (infrared, 
optical) frequencies the inter-band quantum transitions between the electron and hole bands should be taken into 
account and a quantum theory is needed. Recently, a quantum theory of the third-order nonlinear electrodynamic 
response of graphene to a uniform external electric field has been proposed in several publication o^®i®^~®^ds . For 
example, in Refi^ the authors calculated the third-order conductivity a^p^g{uji,uj 2 ,uj 3 ) in the relaxation-free approx¬ 
imation neglecting all scattering effects. In Refi^ the special case of the third-harmonic generation uji = uj 2 = = uj 

has been studied, with the relaxation effects described by a single phenomenological relaxation time. In Refi^ the 
authors generalized their results^ by including the impurities and phonon scattering effects having introduced two 
phenomenological relaxation rates and Te corresponding to the intra- and inter-band electronic transitions. In 
addition, they took into account the finite temperature effects. In the paper— the authors analytically solved the 
semiconductor Bloch equations for the optical graphene response with the electromagnetic interaction described in 
the length gauge {H' = —er ■ E, Ref—). In Ref— they confirmed their analytical results by numerical calculations. 

Although quite a number of theories have been already proposed, their results are still contradictory, apparently 
due to extreme complexity of the problem. For example, the largest contribution to the third harmonic generation 
effect was predicted at fiw ~ 2Ep in Ref— but at Iku ~ 2Ep/’i in Refs— Analyzing the formulas of Ref— we have 
found a number of mistakes (for example, analytical results of Ref— did not reproduce the graphs of that paper) and 
contradictions with Ref—. (After this manuscript has been submitted for publication, the authors of Ref— published 
Frratum, Ref—, in which a misprint in the main formulas of Ref— was corrected.) The development of the theory 
of the third-order nonlinear electrodynamic phenomena in graphene cannot therefore by considered as completed and 
additional independent calculations are required. 

In this paper we develop an analytical theory of the third-order nonlinear effects in graphene and calculate the local 
third-order conductivity To solve the problem we apply a different from Refs— approach, using 

the density matrix formalism and describing the uniform electromagnetic field by a space-dependent scalar potential 
^(r, t) oc exp(iq • r), where the wave-vector q is first assumed to be small but finite, with the limit q 0 being taken 
in final formulas. This way we avoid formally divergent matrix elements of the perturbation with the electron Bloch 
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FIG. 1: Left: The honey-comb lattice of graphene; ai = o(l/2, y/3/2), a .2 = a(—1/2, y/3/2), b = a(0, l/y/l), o is the lattice 
constant. Right: The Brillouin zone of graphene; Gi = (27r/a)(l, l/v^), G 2 = (2Tr/a)(l, —l/\/3), |Gi| = IG 2 I = G = 
47 r/('\/ 3 a); Ki = —K 4 = (27r/a)(l/3,1/V§), K 2 = —K 5 = (27r/a)(2/3, 0), K 3 = —Ke = (27r/a)(l/3, —1/-\/3). Dashed lines 
show the boundaries of the elementary cell in the direct and reciprocal space. 




wave-functions. We study a few specific physical effects, including the third-harmonic generation and the saturable 
absorption effects, and show that the output power emitted by graphene-based nonlinear devices resonantly depends 
on the input frequencies, with the resonance position governed by the electron density in graphene. This allows one to 
control the device operation by the gate voltage which makes the discussed effects very interesting for high-frequency 
electronic and optoelectronic applications. Our results also help to remove the discrepancies between the theories 
proposed so far. 

The paper is organized as follows. In Section |ll] we introduce the tight-binding spectrum and wave functions of 
graphene electrons and calculate matrix elements of some operators which are needed in the rest of the paper. In 
Section Hn] the linear response of graphene^^— is discussed. The results for the first-order conductivity are 

presented in Eqs. - (H51) and in Figure [5] They are used then in Section |VB] 

Section IIVI contains the main results of our work. Here, as well as in Appendixes |B] and o we give a detailed 

/Q\ 

derivation of the third conductivity A summary of results for this function can be found in Eqs. 

(I59II (1681) . In Section|V]we analyze our results in the case of an external monochromatic radiation with the frequency 
(jj and consider two physical effects: the third harmonic generation and the saturable absorption (Kerr effect). Finally, 
in Section EH we summarize all obtained results and discuss the main conclusions from this work. Some technical 
details are collected in Appendixes. 


II. ELECTRONIC SPECTRUM, WAVE FUNCTIONS AND MATRIX ELEMENTS 
A. Tight binding energy and wave fnnctions 

The carbon atoms in graphene occupy a 2D plane (z = 0) and are arranged in a honey-comb lattice, Fig. [Uleft, 
composed out of two triangular sublattices shifted by a vector b with respect to each other. All points of the first 
sublattice (black circles in Fig. [T]a,) are given by the vectors niai -|- n 2 a 2 and those of the second sublattice (open 
circles) by niai -|- n 2 a 2 + b, where ni and n 2 are integers. The basis vectors ai, a 2 and the vector b are chosen as 
shown in Fig. [T^, ai = a(l/2, ^/3/2), a 2 = a(—1/2, •\/3/2), b = a(0,l/-\/3), where a = |ai| = |a 2 | = 2.46 A is the 
lattice constant. 

The carbon atom has four electrons on its outer shell. Three of them, so called cr-electrons, form chemical bonds 
with their neighbors. The fourth (tt-) electron moves in the periodic potential of the honey-comb lattice, 

U{r, z) = Y^ [Ua{r -a,z) + Ua{r - a - b, 2 ;)], 


a 


( 2 ) 
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where Ua is the atomic potential, r = (x, y) is a two-dimensional vector, and the sum is taken over all a vectors of the 
lattice. We will describe its motion within the tight-binding approximation. Within this approach the single-particle 
Hamiltonian is 





( 3 ) 


where p = —ih{dx,dy) and m is the free electron mass. Following the standard procedure of the tight binding 
approximation and assuming only the nearest-neighbors interaction we find the energy^ 

Elk = (-l)'il^kl, (4) 


and the wave functions of the Hamiltonian ([31) 

4';k(r,z) = |;k) = z), (5) 

where 


Mik(r, z) = l^k) = [(-l)'Ck'0a(r - a, z) +'ipair - a - h, z)] 

a 


( 6 ) 


is the Bloch factor, / = 1, 2 is the band index, k = {kx, ky), t is the transfer integral (in graphene t Ri 3 eV), S and A 
are the areas of the sample and of the elementary cell, respectively, and '0a is the atomic wave function. The functions 
Zk and Ck in dH), (H]) are defined as 

Zk = l + = 1 + 2 cos(A:^a/2)e*^'=““/2^ (7) 

Ck = ^k/|^k| = e*'^^ (8) 

where <i)k is the phase of the complex function Z^- They satisfy the relations 

= Zk, Zk+G = Zk, C-k = Ck: Ck+G = Ck, (9) 

where G are the 2D reciprocal lattice vectors which can be chosen as shown in Fig. [T]right, Gi = (27r/a)(l, l/\/3), 
G 2 = (27r/a)(l, —l/v^). The energy Eik and the wave functions |Zk) are periodic in k-space, £'z,k-i-G = Eik and 
10 k -f G) = |/k). The functions (O and (jH]) are normalized as follows 

{l'k'\lk) = [ dr [ dz'l>i,k>{r,z)'liik{r,z) = ( 10 ) 

J S J — CO 


{l'k\lk) = ^ f dr f dzupki'i^, z)uik{r, z) = Sw- 
J A J —OQ 


( 11 ) 


Near the Dirac points, at the corners of the hexagon shaped Brillouin zone k = Kj (j = 1,.., 6 is the valley index. 
Fig. [Ib), the function Z^ assumes the form 


VSa 


7^ Ki -- 

^k- 2 


{-lyH + M 


, kj — k , 


( 12 ) 


where Ki = —K 4 = (27r/a)(l/3, l/v^), K 2 = —K 5 = (27r/a)(2/3,0), K 3 = —Kg = (27r/a)(l/3, —1/-\/3). The energy 
dH) and the function Ck are then (at |kj|a <C 1) 


i^/k = (-l)'^|k,l^(-l)Wft|k,|, 


(13) 


Ckj 


(-l)4^+^fc^ 

|k.l 




(14) 


where the effective (Fermi) velocity is vp = t-\/3a/2/i r; 10® cm/s, and (pk is the polar angle of the vector k, 
k = {kx,ky) = /c(cosi^k,sin(pk). 

In typical graphene structures the chemical potential y, lies in the vicinity of Dirac points, |^| <C t. In this paper 
we will assume that the energy of the photon is also small as compared to the overlap integral t, hio t. Under 
these conditions the energy and the wave-functions of electrons can be described by formulas (fT51) and (imi . Possible 
resonances originated from quantum transitions in other (than Kj) points of the Brillouin zone, as well as the trigonal 
warping effects are not considered in this paper. 
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B. Matrix elements 

To calculate the system response in subsequent chapters we will need some matrix elements with the tight-binding 
functions ©• 


1. Matrix element of the exponential function 


Consider the matrix element (/k|e*'i '’|/'k') of the exponential function e*'! '’. Assuming as above that the neighboring 
wave functions do not overlap we get 


1 


{lk\e^^-'^\l'k') = 4,k'+q2 [1 + (-1) + Ci:Ck-c 


(15) 


The matrix element is nonzero only at k = k' -|- q. Eq. m is valid in the whole Brillouin zone and at arbitrary values 

^Q\ 

of q. Since we aim to calculate the local (q = 0) conductivity cr^^p^g{uji,uj 2 ,uj 3 ) we will need these matrix elements 
at q —>■ 0, however, the terms linear in q must be kept. Then the matrix element (1151) assumes the form 


(;k|e*'i ‘-|^'k')q^o = 4.k'+q (s,i - 


(16) 


where 


dk^ dk^ ■ 


(17) 


Near the Dirac points j = 1 and j = 2 the function Ck is given by Eq. (HI- Substituting this in Eq. m we get 

(18) 


a ■^7’k , i , 

Vkj “ ~^k^ ea/3^/3 


where the upper (lower) sign corresponds to j = 1 (j = 2) and the two-dimensional Levi-Civita symbol is defined in 
m , see Appendix [A] Notice that the more general formula dnl) is valid in the whole Brillouin zone. 


2. Matrix element of the velocity 


Another quantity that we will need below is the diagonal in k (k = k') matrix element of the velocity op¬ 
erator (/k|v|/'k), where v = p/m (here m is the free electron mass). To find it, consider the matrix element 
{lk\[HQ, e^'^ ’^Wl'k'). On the one hand it equals 


On the other hand 


so that 


{lk\[Ho,e^^’^]m = (E,k-Epk')(^k|e*‘i-'-|/'k'). 


[Ho, e*'!-'-] = e*'! '-] = ^ • (ve*'! '- + 


ftq 

T 

In the limit of small q this gives 
hq, 


(/k|ve*'i-- + e^'i-'-virk') = (Aik - A,.k')(^k|e*'i-^|/'k'). 


{lk\{v^,e^^ '^}+\l'k - q) = (Aik - Ei-k-q) f Sm - 


(19) 

( 20 ) 

( 21 ) 

( 22 ) 


where {a, 6}+ = ab + bd is the anticommutator. Now consider two cases, I = V and I ^ V. Eor the diagonal and 
off-diagonal matrix elements we get, keeping only the lowest terms in q: 


qa ( h{lk\va\lk) - 


dEik \ 
dka ) 


= 0 , 


(23) 
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(2h{lk\va,\l'k) - (Ei^ - Ei,k)r^z) =0, I I'. (24) 

Since the direction of the vector q in both formulas (l23l) and (l24t is arbitrary, the terms in parenthesis must vanish. 
Thus we get for the intra- and inter-band matrix elements of the velocity 

(Zk|D„|/k) = 1 = 1 ', (25) 

{lk\v^\l'k) = l^l'. (26) 

The formulas ((25l) and (l26t are valid not only near the Dirac points but in the whole Brillouin zone. Near the Dirac 
point one can use Eq. (fT^ for the energy and Eq. (1181) for the ? 7 ^-function. 


III. LINEAR RESPONSE OF GRAPHENE 


A. Preliminary remarks 

Our ultimate goal is to calculate the current response j(t) of a uniform graphene layer to an external uniform 
electric field E(t) parallel to it. The field is assumed to have an arbitrary time dependence and be presented by 
Fourier expansion 


m = / 

J —( 


E(t) = / dwE^jC 

f —CO 


(27) 


The current response j(t) will have to be calculated in up to the third order in the electric field E. 

The external ac electric field can be included in the Hamiltonian of the system in two ways, using a vector or a 
scalar potential (j){v, t). We choose the second option and write the Hamiltonian in the form 


where the potential can be expanded in the Fourier integral over time, 

/ OO 

-CO 

Notice that at the moment we assume to be space-dependent, 




(28) 


(29) 


(30) 


where the sums are taken over discretized q values, and the transition from the sum to the integral is performed using 
the usual rule 


This has to be done since the field is proportional to the gradient of the potential, E(r, t) = — V(^(r, t), and its Fourier 
component is Eq,^ = —In order to calculate the local (q —0) conductivities and we have, first, to 
find the current response at a finite wave-vector and then take the limit q — >■ 0 in the final formulas. The factor q^iq^j 
should be kept finite in the taking-the-limit procedure. 

In this Section we discuss this process in sufficient detail for the linear response problem. The nonlinear response 
is then treated in Section ITVl 
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B. Density matrix 


The system response to the potential (l29l) is described by Liouville equation for 

where 7 is the phenomenological relaxation rate due to the scattering of electrons 
lattice imperfections (for simplicity, we restrict our consideration by one parameter 
Hq and the equilibrium density matrix po satisfy the eigen-value equations 

H^\\)=E^\X), po|A)=/a|A), (33) 

where |A) = |Zkcr) (we have added the spin variable tr) and f\ = fik is the Fermi distribution function 

fik = (1 + 


the density matrix p 

( 32 ) 

by impurities, phonons and other 
7 ). The unperturbed Hamiltonian 


Expanding the density matrix in powers of the electric field, p = po -I- pi -f p 2 + Ps + • ■ Pn oc E”, and solving the 
Liouville equation in the first order we get 


(A|pi|A')t 



f\' - f\ 

Ex + h{uj + 17 ) 


(A|h^|A')e-*‘"‘, 


(35) 


where = h^{v) = and the subscript t indicates that the density matrix elements (A|pi|A')t depend on 

time [the integrand in (1551) . apart from the exponential function, will then be denoted as (AjpilA')^^]. 


C. Current density 

The two-dimensional current density is calculated according to the formula 

j(ro,t) = -^^(A'|{v5(ro -r)}+|A)(A|p|A')t. 


AA' 


In the first order we get for the a;-Fourier component of the current 




ja;(ro) = Y ^(A'|{v5(ro - r)}+|A) — 


AA' 


Ex' - Ex + h{uj + 17 ) 




Now we expand the potential in the Fourier integral in q, Eq. (1301) . and use the Dirac function representation 


5(ro-r) = 


.iq-(ro-r) 


(36) 


(37) 


(38) 


Then we get 


Jqct. 


e 


E' 


^(A'|{v,e-^U+|A) 


f\' - fx 


AA' 


Ey — Ex + flioj + * 7 ) 


(A|e*'i-'’|A'). 


(39) 


After the substitution |A) = jZkcr) the matrix elements in this expression give the factors (5k',k-q^k,k'-i-q oc 5q_q so that 
the sum over q disappears. 

The formula (1391) is general and valid at arbitrary q. In the limit q ^ 0 we can simplify it, substituting the matrix 
element expansion (I16L thus getting the intra-band (first term) and inter-band (second term) contributions to the 
current: 


Ju) 


ie^gs 


h{uj -f ij)S 


E(^k|ua|Zk) ( - 


Zk 


dfik 

dkp 


ie^gs 

2S 


y^(^k|ua|/k) 


flk ~ flk 


Ik 


Elk - Elk + ^(w + il) 


Vk 


E^. 


(40) 












Here = 2 is the spin degeneracy factor, 


E^. = lim E^. 

q->0 



(41) 


and we have introduced a designation I which means not Z, i.e. 


/2, ifZ = l, 
\ 1, if Z = 2. 


(42) 


Notice that the sums in (l40ll are taken over the whole Brillouin zone of graphene, therefore the valley degeneracy 
factor = 2 has not yet appeared in this expression. 


D. First-order conductivity 

The value in the parenthesis in (1401) is the linear (first-order) local conductivity of graphene 

(1)/ X I* \n\( Y^/fi i» m \ /tk “/ik 

j + -^7—rTT-TT—V'^k- 


h{io + n)S 


/k 


2S 


ik 


- Ei]^ -I- h{u} + ij) 


(43) 


This quantity has been calculated in a number of publications, see e.gi^^— , therefore we do not discuss here further 
details of the integral calculations in (j43|) . Assuming that the temperature is zero, T = 0, and taking into account 
that the main contribution to the sums over k is given by the vicinity of Dirac points we get 




0 ^0.0 ' 


where 


( 1 ) e'^gsgv 

° 16Zi 4Zi 


is the universal optical conductivit} d'^d^ and the dimensionless function 


5, 


(D) = <5„^ (4\Uf^) + 5^ 


(1) 

inter 


consists of two, intra-band and inter-band, contributions: 


(1) 4 i 


TT O -|- iT ’ 


(44) 


(45) 


(46) 


(47) 


41(f^) = zin 


i 2-{n + iT) 


TT 2 -|- (D -|- ir) 


We have introduced here two dimensionless parameters 


0-— r- In. 
H’ H' 


(48) 


(49) 


The intra-band contribution (ITTl) has the Drude form; the inter-band conductivity (H51) has a step-like (logarithmic) 
singularity at huj ~ 2Ep in its real (imaginary) part. At large frequencies hui ^ 2Ep the conductivity assumes the 
universal value (HSl) . Figure [2] shows the intra-band, inter-band and the total conductivity of graphene at F = 0.1. 


IV. THIRD-ORDER RESPONSE OF GRAPHENE 
A. Density matrix and third-order current 

Now we need to calculate the next term in the expansion of the electric current in powers of the electric field. The 
second-order current response oc E^(Z) vanishes due to the central symmetry of graphene. In the third order in the 
perturbation (1^^ the solution of Liouville equation (15^ can be written in the form 
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_,_^^^^_I_^_I_,_^^_ 

0 0.5 1 1.5 2 2.5 3 

Q. 


FIG. 2: (a) The intra-band (I47II . (b) inter-band (I48|l . and (c) total conductivity of a single graphene layer as a function of the 
frequency fl = fajj/Ep at F = h'yjEp = 0.1. In (c) the conductivity in the relaxation-free limit (F = 0, thin curves) is also 
shown. 


/ OO pOO nC 

duji / (ko2 / 

-OO j —OO J — c 




g —z((Vl+tiJ2+<V3)t 


Ex' — Ex -|- h(uji -|- 012 + 
( /a' — /a" 


——- ^ (A|/r.3|A"')(A"'|d.,|A")(A"|V|A') 

wa) -b ihj 


/a" - f\" 


Ex' — Ex'" + h{uJi +bj2) + ih'y I Ex' — Ey + Swi + ifij Ex" — Ey + EjJ2 + ihj 


_ _1_ /_ f\" - /a'" _ f\"' - /a _ 

Ex" — Ex + h{uj 2 -1-013)+ I Ex" — Ex'" + S02 + *^7 Ex'" — Ex + hcus + ih'y 

In order to find the current density in the third-order we should substitute Eq. (1501) in the current definition (1361) . 
Using the delta-function representation (l38|) and the Fourier expansion of the potential (l30l) we get for the third-order 
current: 





doi 


do2 


^4 s 


’f'q2W2 


E 


(A'|{v,e-*+'’}+|A)e*+'’°-d<^i+‘^=+‘^=')* 


E 

A" A" 


qqiq2q3 

(Ale*'!^* I A"') (A"' 16**1= I A") (A" | | A') 1 


AA 


^ Ex> — Ex + h(oi -b 02 + 03) -b ih'j 


f\' - f\" 


f\" - f\" 


Ex' — Ex'" + d(oi 4-02)+ ih'y 


- E 


(Alt 


>|A'")(A" 


m|piq2-r2|;!y//^('y/|p*qiri|x/ 


|A') 


Ex' — Ex" + ftoi -b ih'y Ex" — Ey + fto2 + ih'y 
f\" — fx'" fx'" — fx \ 


A"A" 


Ex" — Ex + h(o2 4-03)4- ih'y 


Ex" — Ex'" + ho2 + ih'y Ey — Ex + hw^ + ih'y 


■ (51) 
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The current density (ISTI) contains the product of three matrix elements of the type (fTSl) . 


(Ale 


,*q 3 rs I ^///^ ^y/l I •r 2 I / X" I ^iqi -ri 


|A")(A" 


|A'). 


(52) 


Each of these factors is the sum of the intra- and inter-band contributions, Eq. (1161) . Expanding the product (15^ we 
obtain a total of eight summands: 

1. One term containing the product of three intraband contributions; we will label it as (3/0) term (three intra- 
and no inter-band factors). This term is purely classical; the corresponding contribution to the current can be 
obtained by solving the classical (Boltzmann) kinetic equation. 

2. Three terms containing the product of two intraband and one interband contributions; we will label them as 
( 2 / 1 ) term (two intra- and one inter-band factors). 

3. Three terms containing the product of one intraband and two interband contributions; they will be labeled as 
( 1 / 2 ) term. 

4. And one term containing the product of three interband contributions; this is a purely quantum contribution; 
it will be labeled as (0/3) term. 

In order to find all these contributions we have to calculate the sums over A = (elk), A' = (e'l'k'). A" = (e"l"k"), 
and A"' = in Eq. (I^Tl) . Perfoming the summations over all spin indexes and all but one factors Ik we 

reduce the expression for the third-order current to the following form (for details see Appendix iBl) : 


( 3 / 0 ) ( 2 / 1 ) ( 1 / 2 ) ( 0 / 3 ) 


(53) 


where 


= 

( 3 / 0 ) 


duJi 


duj2 




u)2 u;3 


• 4 

le gs 


(wi -I- 012 -I- 013 -I- i7)(a;i -I- 012 -|- + ij) h^S 


y^(lk|-0a|lk) 


d^fi 


Ik 


Ik 


dk/sdk^dks 


(54) 


( 2 / 1 ) 


/ du} 

> —00 J —00 


00 nOO 

II duj2 


■'UJl ^UJ2 OJ3 ^ 


E' 


(lk|D„|lk) 


vl 


9 ^ if Ik - fik) 


-ie'^gs 

2S — Eii^ -I- h{uji -I- 012 + W 3 ) -I- ih'y | huii + ih'y h{uJi + 102 ) + ifvy dkpdk^ 


d 


7 

dk 


difik - Ilk) 


hijJi + ih'y dks \ Ej^ — Ei\^ + h{uJi + W2) + ih'y dkp 


d 

dks 


d 


Ei^ — Ei\q -(- h{uj\ -\- C02) T ih'y dk^ ^ ^ik — -1- hwi -t- ihy 


fik - flk 


(55) 


= 

( 1 / 2 ) 


/ CO pCO pOO 

dio, / duJ2 / 

-00 J —00 —00 


ie^gs 






2S h{uji -I- W 2 + W 3 ) + ihy 


vi 


y^(lk|Da|lk) 


Ik 


d 




flk — flk 


d 


Elk — + 1 ^ 2 ) + idy dk^ En^ — Ey^ + huji + ihy ^ 

1 difik-fik)' 


h{uji + UI2) + ihy dks Ei^ — Ejy^ hui -|- ihy ^ 

flk — flk 1 


dZdi 


hwi + ihy Eii^ — Ejy^ -I- h{LUi + UJ2) + ihy dkp 


(56) 
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and 

/ OO pOO pOO 

dio, / dio, / 

-OO j —OO t/ —OO 

( 0 / 3 ) 

ie'^ds _1__ (FklDal^k) _ /3 7 5 

4S' h{uJi + UJ 2 ) + ifij ^ En^ — Ei]^ + h{u!i + W2 + W 3 ) + 

1 1 

Ei^ — i?ik + huji + ih'j Eii^ — Ef^ + fvjJi + ih'y 

One sees a certain structure in formulas (l5l)) - (1571) . The term (3/0) has three poles at low frequencies ujj ~ —* 7 , 
j = 1, 2, 3, and contains the third-order derivative of the distribution function /jt; all other terms have at maximum 
two poles at iOj ~ — 17 . Hence, the term (3/0) is the largest one at low (microwave, terahertz) frequencies. The 
second term ( 2 / 1 ) contains the inter-band energy denominator ^ [Ej]^ — Ei\^ + h(u}\+iij 2 + 0 J^) + ih'f\~"^ and the second 
derivative of the distribution function. Since at low temperatures T/Ep ~ 0 the first derivative 3(/ik ~ flk)/9E 
is proportional to the delta function S{E — Ep), one should expect that the term (2/1) has a second-order pole at 
h(uji +uj 2 + uj^) ~ 2Ep and hence is the largest one at high (mid- and near-infrared) frequencies. The third and fourth 
contributions (1/2) and (0/3) contain the first and zeroth derivatives of the Fermi functions and hence are expected to 
be smaller as compared to the first two terms. It is also worth noting that the (3/0) and (1/2) contributions contain 
the diagonal (intra-band) matrix elements of the velocity (ZklDal/k), while the (2/1) and (0/3) contributions contain 
its off-diagonal (inter-band) matrix elements (/k|{)Q,|Zk). 

Now we define the third-order conductivity and calculate different contributions to it at low (T = 0) temperatures. 



^ ifIV. - fli.)- (57) 


B. Third-order conductivity 


f S') 

The third-order fourth-rank conductivity tensor of graphene W 2 , ^ 3 ) is defined as follows 

/ OO pCO p<X) 

du, / dc.2 / (58) 

-OO J —OO —OO 

From the definition (1551) and the formulas (1571) - ((571) one sees that, first, simultaneous permutations of the arguments 

/ON 

wi, UJ 2 , UJ 3 and the indexes /?, 7 , 6 does not change the third-order current. The conductivity tensor a; 2 , W 3 ) 

can thus be always presented in the symmetric form 

^a/375(wT‘^2,W3) = ^ (wi, W 2 , W 3 ) + (wi, W 3 , Wz) (^2 , Wl, ^ 3 ) 

+ ^i75/3(‘^2,a;3,a;i)-hcfX(‘*^3,wi,a;2)-hdX(w3,W2,wi)), (59) 


where the (unsymmetrical) tilde-conductivities o’^^/^(a;i,a; 2 , W 3 ) can be obtained directly from the expressions (1541) 

- (1571) . Second, in accordance with the expressions (1571) - (1571) the unsymmetrical conductivity W 2 , W 3 ) 

contains four contributions according to the terms (3/0), (2/1), (1/2) and (0/3). We write them in the form 


where 




= (J 


( 3 ) <^( 3 ) 
0 


(Hi, H 2 , H 3 ), 


a 


( 3 ) _ e'^gsgyhv 
° “ 167rF;|, 


2 

F 


(60) 


(61) 


and 


^^ 2 ,1^3) = ^ 2 , n,) + 5X(Hi, H 2 , H 3 ) + + 5X(Hi, H 2 , H 3 ). 


_ <^( 3 / 0 ), 


:(2/i) 


( 1 / 2 ), 


-( 0 / 3 ), 


(62) 


The dimensionless functions (Hi, H 2 , H 3 ) depend on the dimensionless frequencies H^ , j = 1,2,3, and the 

dimensionless scattering parameter F as defined in Eq. (1751) . Introducing the designations 




( 63 ) 
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Oi = 


Hi + zF 


Ol2 = 


Hi “t" H 2 “t“ zF 


O 123 — 


Hi + H 2 + H 3 + zF 


2 2 2 
and calculating the contributions (Hi, H 2 , H 3 ) (see Appendix [Cl for details) we get: 


c(3/o)/q q O';— »^op7o 

80^^^0i20l’ 

I ^aS^/S'y “t“ (^0:^75/4)0123 ^a^'yS/ 

4 y OiOi2(l + 0123)^ Oi(l + Oi2)(l + O123) Oi(l + Oi2)(l + 0123)^ 

— (SajsSjS — Aq^^ 5 / 4 )J 7 i(Oi, O12, O123) + {SajSj 3 S — J^2(Oi, O12, O123) ) 


iAoi^^s 




(64) 


(65) 


— O 12 , 0123 ) —>■ (—Oi, —O 12 , —O 123 )), 


( 66 ) 


5i'/7i(^i’^2,H3) = j 


/I _|_ ^ '\ ^ 0 / 375/4 (‘^a^^75 ^ 0 / 375 / 4 ) 

O123O12O1 0 l 230 l(l + O12) 


4 y O 1 O 123 O 12 

^a ^'^75 ^ 0 / 375 /4 


J 3 (Oi, O12) + 


*^ 07^/35 ^^ 0 ^ 75/4 


Oi23 ^ ^ O 123 

~((Oi, O12,0123) —>■ (—Oi, —O12, —O123)), 


Ji{Oi,Oi2) 


‘^i°/7i(^l’^2,H3) = -r-^^j5(Ol,Ol23)- 




320 


12 


(67) 


( 68 ) 


In the third lines of Eqs. (1661) and (1671) one should subtract the corresponding expressions from the first two lines, 
with the arguments (Oi, O 12 , O 123 ) being replaced by (—Oi, — O 12 , — Oi 23 )- The integrals Jn in (l65|) - (l 68 l) are defined 
as follows 


Ji{a,b,c) = / 


dx 


1 1 
+ 


1 


1 (x + a)(a; + 6)(a; + c) \ a;^ x(x + c) {x + b)(x + c) (x + c)^!’ 


(69) 


a;(a; + a)(x + b)^(x + c)'' 
°° dx 


/ OO 

J73(u, 6) / 2 / \/ I Li ’ 

7 i x‘‘{x + a)[x + b) 

^ x{x + a)(a: + 6)2 ’ 

dx f 1 1 \ / 1 


^ dx f 


X + a X — a I \ X + b x — b 


Their explicit expressions are: 


, , 1 —ab + 2ac + 36c — 4c^ 

Ji(a, 6,c) = — + 


abc c(c — a)2(6 — c)^(c + 1) (a — 6)(6 — c)2(6 + 1) (a — c)(6 — c)(c + 1)^ 

—— 2a^c + 3a6c + ac^ — 6c^ —2a6 + ac + 36^ — 6c 

2 ? 712 z 13 ln(u T 1) j^ 2 z u\2{u 72 m(6 -t- 1) 

+ a6 — 4ac — 36c + 5c^ 


b‘^(a — hY{b — c)2 


c(a — c)3(c — 6)2 


ln(c + 1 ), 


(70) 

(71) 

(72) 

(73) 


(74) 
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J2{a,b,c) = 


1 


b{a — b){b — c){b + 1 ) a{a — b)^{a — c) 


ln(a + 1 ) 


2 ab — ac— 36^ + 2 bc 
b'^{a — byih — cY 


ln (6 + 1 ) - 


1 


c(c — a) {b — cY 


ln(c + 1 ), 


JYa,b) = ^ + —— ln(l + a) - —— ln(l + 6 ), 


ab a?-(a — b) 


JA{a,b) = - 


b{a — b){l + b) a(a — b) 


b^(a — b) 


ln(l + a) + - —T ln(l + b), 


b^{a-bY 


(75) 

(76) 

(77) 


775(0,6) 


4 2 b j (1-a) 

ab a?{a — h){a-\-b) (1 + a) 


2 a 1 (1 - 

62(a-6)(a + 6) 


(78) 


The formulas (1741) - (|78l) are valid at arbitrary values of the arguments a, b and c, but are not always convenient for 
numerical plotting of the results since the factors (a —6), (o —c), (b — c) in the denominators may cause error messages 
“division by zero”. It is easy to verify that in all limits a ^ b, a ^ c, b ^ c, etc., the integrals Jn remain finite [this 
is seen directly from the definitions (l69ll - (l7^ : notice also that the arguments a, b, c are complex at a finite F]. The 
explicit expressions for Jn in the special cases a = b, a = c, etc., are provided in Appendix IdI 

As we have discussed above, the contribution (3/0) has three poles at low frequencies [the factors I/O 123 O 12 O 1 in 
Eq. (l65]) ]. and the contribution (2/1) has the second-order pole at h{uJi +UJ 2 +UJ 3 ) — ±2Ep [the factors 1/(1 ± 0123 )^ 
in Eq. (1551) and 1/(1 ± cY in Eq. (1751) ]. In addition, the function has the first-order pole at the two-photon 

absorption frequency h(uji +UJ 2 ) — ±2Ef [the factors 1/(1 ± O 12 ) in Eqs. (1551) and (1571) and 1/(1 ± b) in Eqs. (175)) - 
(EH], as well as (weaker) logarithmic singularities at O 123 ~ ±1, O 12 ~ ±1, and Oi ~ ±1. 

The third-order conductivity tensor (l59l) satisfies the relation 



= a 


(3) 

a/S'yS 


(W1,W2, W3). 


It has three independent components 


= ct(3) 

^ xyxy ^yxyx'> 


= ct(3) 

^ xyyx ^yxxy^ 


= ct(3) 

^ xxyy ^yyxxi 


and the component cr^Jxx = CTyyyy which satisfies the equality 


a(3) = 

xxxx ^xxyy 


.(3) 


7 ( 3 ) 

xyxy 


(79) 


(80) 


(81) 


This agrees with the phenomenological relations for the third-order response function in an isotropic medium^. 

The formulas obtained in this Section are the main result of this work. They show that the Fermi-energy related 
resonance in the third-order conductivity is the case at fli -I- 172 + II 3 ~ 2 (which corresponds to huj = 2 Ef/2> if all 
three frequency arguments are the same), as it was first pointed out in Ref.—. In our previous publicatioit^^ not 
all contributions to the mixed current terms (2/1), Eq. (1551) . and (1/2), Eq. (I56L have been taken into account; 
as a consequence, the main resonance in the third-harmonic-generation term was found at huj = 2 Ef instead of at 
fvjj = 2EF/'i and the absolute value of the third-order conductivity was overestimated. In the relaxation-free limit 
r —>■ 0 the results (155)) - (I68|) mainly agree with those obtained in Ref—, with the exception of a small vicinity of the 
pole at [ 171 - 1 - 172 - 1-173 — 21 < F where a certain discussion is needed (it can be found below in Section IV Al) . The 
results (I5ni)-(I551), which have been obtained by the method different from that used in Ref—, also agree with those 
of Ref.— if to correct them according to the Erratum^. 

In the rest of the paper we discuss some consequences from the obtained results. In particular, the corrected results 
for the third-harmonic generation effect are considered in Section IV Al 


C. Asymptotes at low and high frequencies 

First, consider the system response at low frequencies, when all frequency parameters satisfy the conditions |I7j -|- 
iF| <C 1, or rij ^ 1 and F = h/EpT <C 1. The condition F 1 means ypT^im^ ^ 1, i.e., the mean-free path 
should exceed the average distance between electrons. At r = (0.1 — 1) ps it means Ug > 10^° cm“^. Since a typical 
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FIG. 3: The relations hio — 2 Ef = 2hvFy/En2 (label 2) and huj = 2Ef/3 (label 2/3) plotted as the frequency / = a;/27r and 
the wavelength A = c/f versus the electron density Us- 


graphene electron density exceeds ^ 10^^ cm“^, the condition F <C 1 is usually satisfied in real experiments. The 
second condition, hio <C Ep, is satisfied at up to mid-infrared frequencies, see Figure |31 For example, at the density 
ns ~ 10^^ cm“^ this condition corresponds to / <C 30 THz. 

/ox 

Under the conditions Iflj -|-ir| <C 1, j = 1,2,3, the term (3/0) gives the largest contribution to a^p^g{uJi,u} 2 ,uJ 3 ), 
and 


(Wl, W2, W3 ) ~ (Wi, W2 , W3 ) = 


(3) 


.(3) 


iA 




6{ni + iT){n2 + iT){n3 + ir)' 


-|- ir| <c 1. 


(82) 


( 3 / 0 ) 


This result was obtained within the quasi-classical approach in Ref.—, and then within the full quantum theory in 
Refs.— It describes the low-frequency (microwave, terahertz) response of doped graphene. 

In the opposite limit, when all frequencies are large, IIj ^ 1, the most important contributions to cu 2 , W 3 ) 

come from the imaginary parts of logarithm functions in (IMl) - (l 68 ll . (IT^ - (I78L 


lim In 

r 2 i—)-CO 


1-bOi 

l-Oi 


= in 


(83) 


(and similarly for O 12 and O 123 ). This gives 

The result 


ttA, 




r< 1 < |Oj|. 


1 , 2 , 0 3(0,+f}2)(02 + F! 3 )(f 23 + 0i)(F!i+02 + f 23 )’ 

agrees with the one obtained in Ref—. It describes the optical response of intrinsic graphene. 


(84) 


V. THIRD-ORDER RESPONSE TO A MONOCHROMATIC RADIATION 


Now we consider the third-order response of an isolated, freely hanging in air graphene layer to an external monochro¬ 
matic radiation with the frequency uj and the intensity Iu> normally incident on the layer. Let the external electro¬ 
magnetic radiation be linearly polarized in the a;-direction and have only one frequency harmonic, 

Ep{t) = EqSxp cosijt. (85) 


Then the Fourier component of the field is 

~ '^EqSxP ^^(wi — w) -|- i5(a;i -I- tu)^, 
and the intensity of the incident wave reads 




( 86 ) 


( 87 ) 
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FIG. 4: Four contributions to the third order dimensionless conductivity (a) (3/0) contribution, (b) (2/1) 

contribution, (c) (1/2) contribution, (d) (0/3) contribution. The effective relaxation rate is F = 0.1. 


Eui = Eo/2. The electric field really acting on the graphene electrons Eo/{l + 27rcr(^)/c)] differs from (IM|) due to 
the self-consistent screening effects (for their influence on the third-order response of graphene see, e.g., Refi^). The 
difference is determined by the parameter 2'k<7^^'^ (w)/c which is small (< (1 — 2)%) at infrared and optical frequencies. 
If to ignore this difference the third-order current (1581) assumes the form 

e-“*-hc.c.^, (88) 

/'Q'j 

where c.c. means the complex conjugate. Since ayJxx equals zero, the current is polarized in the same {x-) direction 
and has two contributions, one at the third harmonic Sw and another at the frequency of the incident wave w (the 
Kerr effect). 




1^0 


r(3) 


-(- (cu, CU, (jj^xxxi^^ ^ T E (jj^xxxi 


A. Third harmonic generation 

The third-harmonic current in Eq. (1881) with the Fourier component = ai^xx(uJ,uj,uj)(Eo/2)^ is determined by 
the conductivity tTxxxx(w, w, w). The Fourier components of the induced electric and magnetic fields at the frequency 
3w are 


Esc^ = = ^crEL(w,w,w)F^. 

The intensity of the emitted third-harmonic wave (in the positive and negative z-directions) is then 


(89) 


hu, = = I — 


^XXXX (^5 OJ, Oj) 


( 90 ) 
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FIG. 5: The total third order conductivity under the conditions of Figure 3(b) of Refi^. The parameter 

F = 0.11 corresponds to the numbers (the relaxation rate 33 meV, the Fermi energy 0.3 eV) used in Ref.—. 


it can be written as 


77 


\h'^c^ Trn^hvp 




(91) 


where Si^Jxxi^, is a sum of four contributions, Eq. (15^ . These contributions are shown in Figure|3](at F = 0.1). 
The following features are seen: 

1. The classical (3/0) contribution has a large resonance at low frequencies fl < F. The contributions (2/1) has 
a sharp resonance at fl ~ 2/3 corresponding to the three-photon absorption edge and weaker singularities at 
n ~ F and at SI ~ 1; the latter corresponds to the two-photon absorption. The contribution (1/2) has resonances 
at SI < F, SI ~ 1 and SI ~ 2, and the purely quantum contribution (0/3) has weak (logarithmic) singularities at 
the frequencies of the single-photon (SI ~ 2) and triple-photon (SI ~ 2/3) inter-band transitions. 

2. The absolute values of the different contributions to the third-order conductivity gradually decrease from the 
purely classical (3/0) to the purely quantum (0/3) contributions: at F = 0.1 the scale of the (3/0) term is about 
few thousands, of the (2/1) term - about 10, of the (1/2) term about 1 and of the (0/3) contribution ~ about 
0.1. The largest contribution at low (microwave, terahertz) frequencies is thus expected from the classical term 
(3/0), while the largest (and resonant) contribution at higher (mid- and near-infrared) frequencies is expected 
at the resonance SI = 2/3 {hui = 2Ep/3) from the (2/1) contribution. 

Figure E] shows the total third-harmonic conductivity 5i;7xx(Sl, S2, SI) calculated with our formulas (l54ll - (|57)) at pa¬ 
rameters corresponding to Figure 3(b) from Refi^. The figures agree with each other. 

The most interesting feature in Figures 0] and [S] is the case near the triple-photon absorption edge at SI = 2/3 (it 
was not discussed in detail in Ref.—). The behavior of the third-order conductivity near this point, at small values of 
F and especially in the limit F —0, is rather nontrivial. To illustrate this we show in Figure [5] the total dimensionless 

/q\ 

third-order conductivity 5/xxx(Sl, SI, SI) (the sum of all four contributions) near the point SI = 2/3 and, for comparison, 
near the point SI = 1. The real and imaginary parts of this function are shown at F = 0 (the relaxation-free limit, 
black curves) and at F = 0.01 (red curves). Consider, first, the vicinity of the point SI = 1 (the two-photon absorption 
edge). Figure EDb). If F = 0, the real part of Sx^Jxxi^,^,^) has a step-like feature near the point SI = 1, while the 
imaginary part - a logarithmic singularity, similar to those which are the case in the linear response conductivity. Fig. 
[2j At the finite value of F = 0.01 these singularities are smeared out, as one would expect. Near the point SI = 2/3 

f3l 

the behavior of the real and imaginary parts of 5ixa:x(Sl, SI, SI) is completely different. At F = 0 one sees, again, 
the step-like and logarithmic singularities in the real and imaginary parts, respectively. At F = 0.01, however, these 
singularities are not smeared out but an additional, much stronger, resonance appears near SI = 2/3. To understand 
this behavior we have analyzed the general formulas (15^ - (1551) and found that near the point SI = Slg = 2/3, apart 
from the logarithmic singularity, the function Sxxxx{^,^,^) has the contribution 


3 F 


^x^xxxi^j 


32 (Sl-Slo + ir/3)2’ 


at |Sl-2/3| <r. 


(92) 
















FIG. 6: The total third order conductivity (a) near the triple-photon absorption edge 3hfl ~ 2Ep, and (b) 

near the double-photon absorption edge 2hQ ~ 2Ep, in the relaxation-free limit F = 0 and at a finite F = 0.01. In the panels 
(c) and (d) the exact solution (I62II is compared to the approximation (1921) at (c) F = 0.005 and F = 0.001. The width of the 
resonance decreases, its height increases, and the accuracy of the approximation ((9^ improves when F tends to zero. 


which results from the second order poles ~ (l±Oi 23 ) ^ in [and ~ (l±c) ^ in (174)) ]. How well this approximation 
reproduces the exact expression for is illustrated in Figures [ 6 l[c),(d). One sees that the smaller F, the 

closer is the approximation ((Ml) to the exact expression, and the larger is the resonance. When F decreases, the width 
of the resonance tends to zero while its height tends to infinity since exactly at 0 = 2/3 Sxxxx{^, O, O) oc 1/F. 

In Refi^ the authors argued that their results (obtained at finite values of the relaxation parameters) recover those 
of Reff^ in the limit Fg, F^ —>■ 0. This is however valid only at frequencies far from the resonance (EH). As seen from 
the above discussion the resonant contribution to Sx^Jxxi^^,^^,^), Eq. (I92L could not be obtained in the relaxation- 
free theory^^, since the relaxation parameter F was set to zero from the very beginning. It should be emphasized 
that, since in real experimental systems the F parameter is small (F <C 1) but always finite, the contribution (1921) is 
dominant around the frequency fl ~ 2/3. The position of this resonance depends on the Fermi energy, therefore this 
resonance can be used for creation of efficient, resonant, voltage tunable optoelectronic devices operating in the mid- 
and near-infrared spectral range. 

Now consider the up-conversion efficiency of a single graphene layer at low and high frequencies. In Figure 0 
we show the intensity of the third harmonic, normalized to the intensity of the incident wave Eq. (1911) . as 
a function of the electron density. At infrared frequencies. Figure [71[a), the up-conversion efficiency is of order of 
10“^ — 10“®, dependent on the wavelength, at the chosen set of parameters {luj = 10 MW/cm^, /17 = 1 meV) and 
resonantly increases around the points huj/Ep = 2, 1, and 2/3. These resonances are very narrow and sharp which 
can be used for a fine electrical tuning of the graphene based frequency multipliers. At terahertz frequencies. Figure 
Elb), the up-conversion efficiency is substantially higher (for this graph we have chosen a ten times lower input power 
density = 1 MW/cm^ and a ten times larger relaxation rate Ivy = 10 meV) and smoothly depends on the electron 
density. 
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FIG. 7: The up-conversion efficiency Isui/Iw, Eq. (I91II . as a function of the electron density, in the range of (a) infrared and (b) 
terahertz frequencies. A is the wavelength of the incident wave. The resonances in (a) correspond to 12 = 2, 1, and 2/3 (from 
left to right). 


The results obtained in this Section are valid for an isolated graphene layer freely hanging in air. If graphene lies on 
a dielectric substrate metalized from the back side the intensity of the third harmonic emitted from such a structure 
can be more than two orders of magnitude larger, see Refi^. 


B. Absorption saturation 


The third-order response of graphene at the frequency of the incident wave w is given by the terms in Eq. 
proportional to Using the symmetry of (t^Jxx we get 




73 ^( 3 ) 


H/Q U XXXX 


■uj, u), u) e 


+ C.C.. 


(93) 


Figure [5] shows four contributions to the third order conductivity Six^xx(—^, fl). One sees that, at low frequencies, 
n < r, the largest contribution is, again, the classical one, (3/0). It is of the same order of magnitude as the 
function six^xx(^,^,^), compare with FigurelH At higher frequencies, 17 > 1, the largest values of iSxxxx(—17,17,17) 
are achieved at the condition 17 Ri 2 which corresponds, again, to the second-order pole at 17 1 -|-172 -I- 173 = 2. The 
(2/1), (1/2), and (0/3) contributions are of the same order of magnitude and are much larger than the corresponding 
contributions to 5Ma:a:(17,17,17). The term s'^Jxl (—17,17,17) has the second-order pole at 17 ~ 2, see Figure [5Kb); the 
(1/2) and (0/3) terms have logarithmic singularities near this point. 

Figure UKa) shows the total third-order dimensionless conductivity 17,17,17) as a function of the frequency 

17 under the conditions of Figure 4(b) from Refj^ (at F = 0.11 which corresponds to the relaxation rate 33 meV and 
the Fermi energy 0.3 eV used in Ref<^). The figures agree with each other. In Figure IHKb) we show the vicinity of 
the most important feature at 17 ~ 2 and at T = 0.01. The frequency dependence of the real and imaginary parts of 
17,17,17) is quite unusual. This is a consequence of the crossing and interaction of three resonant singularities, 
a second order pole and two logarithmic singularities. This is illustrated in Figures|nKb,c). In these graphs we show the 
real (Figure IHKb)) and imaginary (Figure IHKc)) parts of the function 17i, 17,17) where the first argument —17i 

slightly differs from —17. Consider, for example, the black curves in Figures |9Kb,c) which correspond to 17i = 0.917. 
One sees that these curves have two weak logarithmic singularities (a step-like behavior in the real part and a weak 
resonance in the imaginary part) at the frequencies 17 = 2 and 17 k, 2.22 (the latter corresponds to 17i = 0.917 = 2). 
In addition, they have a strong second-order pole at 17 Ri 1.82, which corresponds to the three-photon absorption 
frequency 217 — 17i = 2, i.e., 17 = 2/1.1 k, 1.82. All these resonances have a small amplitude (< 40). When 17i 
approaches 17, see the red dotted curves for 17i = 0.9517, the pole is shifted to the point 17 = 2/1.05 r; 1.9 and the 
logarithmic singularity 17i = 2 to the point 17 = 2/0.95 Ri 2.1. Since all three resonances get closer their amplitude 
increases. If 17i = 17, Figure (HKa), all three singularities meet in one point 17 = 2 and three small resonances (with 
the amplitudes < 100) create a large resonant structure with the amplitude ~ 2000 — 4000. When 17i becomes larger 
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FIG. 8: Four contributions to the third order conductivity (a) (3/0) contribution, (b) (2/1) contribution, 

(c) (1/2) contribution, (d) (0/3) contribution. The effective relaxation rate is F = 0.1. 



than (= 1.05, 1.1; green and blue curves in Figure inKb,c)) the three resonances move away from each other and 
their amplitudes are reduced again. 

The third-order current contribution (IMl) can be added to the linear-response current, Eq. (|4^ . to get the electric- 
field dependent response of graphene at the frequency w. In up to the third order in the electric field we have 


jx,uj(t) ~ ^Eq 




e + c.c. =-aooxii^,Eo)Eoe “*-|-c.c.. 


(94) 


The corresponding field-dependent dynamic conductivity is then 


axx{uJ,Eo) = crW(w) -h-E^cr®^^(-w,w,w) 


.( 1 ) 




(95) 


where the functions 5^^^(f2) and fl 2 , fis) are defined in (H51) - (HSl) and (15^ . (1551 - (IB51) . respectively. The 

dimensionless field parameter E = cEq/E pkp in (1951) is the work of the field, performed on 2D electrons at the length 
~ kp^, divided by their average (Fermi) energy. Figure ITOl shows the absorption coefficient 


A{uj,Eo) 


|1 -I- 2TTaxx{uJ,Eo)/c\'^ 


(96) 


of an isolated graphene layer with the field-dependent conductivity (1551) at low (D ^ 1) and high (D ~ 2) frequencies 
(cr' is the real part of cr; notice that the self-consistent-field effect is taken into account here). The absorption is seen 
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FIG. 9: The third-order dimensionless conductivity as a function of the frequency Q, (a) under the conditions 

of Figure 4(b) from Ref.— (F = 0.11) and (b) in the vicinity of the most important feature at 11 ~ 2 and F = 0.01. The real 
(c) and imaginary (d) parts of the function Hi, H) for Hi = 0.9H, 0.95H, 1.05H and I.IH. One sees that the large 

(amplitude ~ 4000) resonant singularity in (a) is a consequence of the crossing and interaction of three smaller (amplitudes 
~ 100) resonances: the second-order pole at 2H — Hi = 2, and two logarithmic singularities at H = 2 and Hi = 2. 


to be substantially suppressed at the external electric field parameter exceeding a certain value which depends on 
the scattering rate F. If F ~ 0.1, the absorption is noticeably suppressed at F < 0.05 — 0.1. At the electron density 
ns ~ 10 ^^ cm“^ this corresponds to the electric fields 

Eo ~ (0.05 - 0 ^ (10 - 20) kV/cm. (97) 

e 

If F is ten times smaller, Figure fTOf c). the required electric field is reduced by, roughly, one order of magnitude. It 
is worth noting that at small values of F the resonance near H ~ 2 is very narrow which can be used for a precise 
control of the saturable absorption of infrared-laser mirrors, cf—i^. 


VI. SUMMARY AND CONCLUSIONS 

We have presented a quantum theory of the third-order nonlinear electrodynamic response of graphene. The 
obtained results can be used for analysis of a large number of different nonlinear phenomena, such as the harmonics 
generation, the four wave mixing, the current induced sum and difference generation, and other. 

The results show that there are two main frequency ranges where the nonlinear effects in graphene are especially 
large and therefore interesting for applications. The first range corresponds to low (microwave, terahertz) frequencies 
described by the condition H < F. The second range corresponds to the three-photon absorption-edge resonance at 
high (infrared) frequencies (Hi -|- H 2 -I- H 3 — 2) < F. There also exist two-photon and one-photon absorption-edge 
singularities which are, however, much weaker. The positions of high-frequency resonances depend on the electron 
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£1 

FIG. 10: The dimensionless absorption coefficient (1961) in a single (isolated) graphene layer at (a) low and (b) high frequencies, 
r = 0.1, and (c) high frequencies at F = 0.01. The curves are labeled by the value of the dimensionless electric field 
F = eEo/EpkF. 


density and can therefore be electrically tuned in a broad frequency range, thus opening interesting opportunities for 
electrically tunable device applications. 

The meaning of the relaxation parameter F may be different in these frequency ranges. At low frequencies T 
corresponds to the intra-band (Drude) scattering rate, F fijr. In currently available graphene samples the scattering 
time T is about 0.1 — 1 ps, so that at low frequencies the strong nonlinear phenomena in graphene are to be expected 
in the technologically important range of microwave and terahertz frequencies f ^ W THz. At infrared frequencies 
the width of the three-photon absorption resonance is determined by the inter-band relaxation parameter F, which 
may differ from the intra-band one. We do not analyze the resonance linewidths further in this work since it would 
go far beyond its scope. The F parameters introduced here should just be treated as phenomenological quantities 
to be extracted from the experiment. Their dependence on the electron density, temperature, excitation power (the 
electron heating effects) may be investigated in future when enough experimental data will be accumulated. 

Our work removes contradictions between results of different theories which have been published so far. This work 
corrects our results published in Refi^ and is now in agreement with the relaxation-free theory of Refi^ (with the 
reservation concerning the second-order-pole contribution (IMl) discussed in Section IV Al) and with the full theory of 
Refi^ whose analytical formulas have been recently corrected in Refi^. 

A comparison with experimental data is quite difficult at the present stage. There are several reasons for that. 
First, most of the experiments have been done at a single (fixed) frequency or in a narrow frequency range, often 
at the optical frequencies. The value of the Fermi energy and the effective scattering rate 7 is usually unknown in 
such publications, which impedes the quantitative comparison with the theory. Second, in experiments graphene 
usually lies on a substrate and is sometimes covered by another dielectric material, while the theories treat a single 
graphene layer freely hanging in air. It was shown in Refi^ that even in a simple structure with graphene lying on 
a dielectric substrate with the metalized back side the third-order response can be both several order of magnitude 
larger and several orders of magnitude smaller than in the isolated graphene layer. A precise knowledge of graphene 
and substrate properties is thus needed for a quantitative comparison of the theory and experiments. Third, many 
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nonlinear experiments are made in the pulsed-excitation regime, with the power density level ^GW/cm^, which lies 
far beyond the applicability area of the perturbation theory used in our work and in other paper a^^i^^ . The nonlinear 
electrodynamics of graphene is thus still at the initial stage, and further experimental and theoretical studies of this 
interesting area of the fundamental and applied physics is highly desirable. 
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Appendix A: Some useful formulas 

Some definitions and formulas used in the main text are presented here. 
The two-dimensional Levi-Civita symbol is defined as 


— 


0 1 
-1 0 y ■ 


(Al) 


It satisfies the following relations useful by calculating the integrals in Appendix O 


ya + SxJ 

fl(T H” ^Xa^fiy^ — ^a/3^7(5 “1“ ^ a-j^/SS H” ^aS^l3'y- 

The integration over the polar angle in the k-space is performed with the help of the formulas 

U2 


1 

(J^a.ki^') = ^ J kctk^d(l) — 


(A2) 


(A3) 


1 / \ 

{ka.k^kjkg^ = J kQ^k^k^k^d(l) — ^^0/3^75 H" do^S^/S'yj ■ 


(A4) 


Appendix B: Derivation of formulas (1531) — (I57D 


To simplify the current expression (1511) we first perform the following permutations of the integration variables. In 
the last line of Eq. (EU we replace (q 2 W 2 r 2 ) O (qiwiri) and then (qswara) (q 2 a; 2 r 2 )- This gives: 


i(3) 


/ oo roo POO ^4 

dwi / duj2 / 

-OO ^—oc */— oo _ 


E 


(A'|{v, 


qqiq2q3 aa' 


E 

A"A" 


- E 


A"A" 


- E 


A"A" 


+ E 


A"A" 


(A|e^q^ fx- - fx" 

Ex' — Ex'" + h{oJi J- a;2) J- Ex' — Ex" -I- fujJi + ih'y 

(A|e^q=^'^=^|A"')(A'"|e^q^-^^|A")(A"|e^qi'^nA') fx" - fx'" 

Ex' — Ex'" + fi{uji + UJ 2 ) + itij Ex" — Ex'" + fiw2 + *^7 

(A|e^q^'^^|A"')(A"'|e^qi‘•HA")(A"|e^q^--^lA') fx" - fx'" 

Ex" — Ex + h{uji +UJ2) + idj Ex" — Ex'" J- fvjJi + ihj 

(A|e^q^'^^|A'")(A"'|e^qi‘•HA")(A"|e^q^ ‘-^|A') fx'" - fx 

Ex" — Ex + fi{uji +UJ 2 ) + ih'y Ex'" — Ex + huj2 + *^7 


(Bl) 
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Then in the second and forth lines in the square parenthesis we replace (q 2 W 2 r 2 ) (qiwiri): 

„4 (A'|{v, 


j(^^(ro,t) = / duji duj2 / dujs— ^ 


' —CO «/ —CO 


E' 


qqiq2q3 aa' 


E 

A"A" 


(A|e^q^‘•^|A"')(A"'|e^q^'^^|A")(A"|e^q^'"^|A') /a> - fy 

Ey — Ex'" + + L 02 ) + ifij Ex' — Ex" + fnoi + ihj 

(A|e^q=’'^=’|A" 0 (A"'|e^qi'^H-^")(A"|e^q=^-^^|A') fx" - fx’" 

Ey - Ey" + h{uJi + W 2 ) + ih'y Ex" - Ey + ftwi + ih'j 

(A|e*‘i= '-=|A"')(A"'|e*'ii ‘-nA")(A"|e*'i3 ‘-3|A') fy - fy 


- E 

A"A" 

- E 

A"A" 


E 

A"A" 


Ey — Ex + fi{uJi + UJ 2 ) + Ey — Ey + fiwi + ifiq 

(A|e^q^'^nA"')(A"'|e^q^ ‘'^|A")(A"|e^q^ ‘'^|A') fy - fx 

Ey — Ex + fi{uji + UJ 2 ) + ih'y Ex'" — Ex + huji + ih-f 


(B2) 


Now in all denominators we have the same frequency factors, wi + 17 , wi + a ;2 + *7 and wi + u ;2 + wa + 17 . Substituting 
the matrix elements (fTOl) at small q and taking the sums over the spin indexes and over the momenta k'", k", k' (we 
use the momentum conservation factors dk".k'+qi, etc.) we get 


/ OO pOO p 

duJi / duj2 / 
-00 J —00 J — 


°° . 


i: 


i-’qicji H^c[2^2 V^qa^^s ’• 


pi(q3+qi+q2)-ro-'i(a;i+cj2+^^3)* 


E 


(/'k|{v, e-'(q3+qi+q2) r|^|;^ k + qa + qi + q 2 ) 


7^ El']ii — i?i^k+q3+qi+q2 + ^(^1 + 022 + Wa) + *^7 


WVi 


E 


Ow" - 


2 935^k+q2+qi 


(- 1 )' 


l"+l' 


<727?7k+qi ) ( Si"I' - ^ ' qi/3Vu 


El'lfi — £’i'"^k+qi+q2 + E{lU\ + UI 2 ) + ih'y 


fl'k — /i",k+qi 


El'k .^/",k+qi T E 01 -f ih'j 


- E 


SW" - -g35?7k+q2+qi ) ( Sl'"l 


(-if 


gi/3?7k+q, ) ( Sl"l' - ^ ^^2 ' g27^k 


l"+l' 




El'k — i?i'",k+q2+qi + h{uJi + UJ 2 ) + ih'j 

/i",k+q 2 ~ /i"',k+q 2 +qi 


.^/“,k+q 2 .^/" 7 k+q 2 +qi huJi -f ih'j 


- E 


Sw" - ^ g27??k+q3+qi ) (Si'-i" - ^ '^‘2 di0Vk+q,) (si'n' - ^ 




-Q36Vw 


l"l'‘ 


£'i",k+q3 ~ £’z,k+q3+qi+q2 + h{uJi + UJ 2 ) + ih'j 


/i",k+q3 /r",k+q3+qi 


.^/“,k+q 3 .^/" 7 k+q 3 +qi T hiO\ -f ih'j 


E 

i"i"' 


Sw" - ^^^'7l/3?if+q3+qJ lSl'"l" - ^^4^927^k+q3 ) ^ 




-Q36Vw 


El" .^/,k+q3+qi+q2 ~b ^(cCl -f CU 2 ) “t“ ih'j 


_/i"',k+q3+q2 /z,k+q3+qi+q2_ 

k+q3+q2 ~ £'/,k+q3+qi+q2 + hlOl + ih'j 


(B3) 


Now consider different contributions to the current. 
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1. The contribution (3/0) 


The purely intra-band contribution (3/0) is obtained from Eq. (IB3I) if in all matrix elements only the Kronecker- 
symbol terms taken into account. Taking the sums over Z'", and I' we get 


(3/0) 


du!i 


f 


dujo 


duj; 


4 

e ffs 

2S 


E 


^qiLOi 'rq2<^2 'rqa^^s 


^i(q3+qi+q2)-ro-i((^i+(^2+^3)i 


E 

Zk 


qiq2q3 

{lk\{va, k + qa -f qi -f q2) 

Eik — E';,k+q 3 -i-qi-(-q 2 + Zi(wi -f a ;2 “b wa) + itvy 


fl\^ /"tk+qi 


Elk - £i/,k+qi-i-q 2 + + 0 J 2 ) + ih'y Elk - Ei,k+q^ +Huji+ ih'y 

_ __ fl.k+q2 ~ /tk-tq2+qi _ 

Elk — El^k+q 2 +m E h{uJl + UJ 2 ) + itv^ El,k+ci 2 ~ ^i,k+q 2 +qi + ^1 + 

_ __ /tk+qa ~ /i,k+q3+qi _ 

El,k-\-ci3 T^tk+q3+qi+q2 T h,{u)\ “b ^^ 2 ) “t“ /Zty E/ k+qs .ZD/ k+qa+qi “t“ “b iHrf 

_ /i,k+q3+q2 ~ /tk+qa-Tqi+qa _ 

ZD/.k+qa ZD/ k+qa+qi+qa T Zi(cCa “b OJ 2 ) ~b ih'J iD/ k+qa-t-qa ZD/ k+qa+qi+qa “b Eol -b ih'y 
One sees that all the Fermi-function differences are proportional to qi and can be presented as 

dfik 


fik - fi. 


k+qi 




dkp ’ 


(B4) 


(B5) 


and similarly for the second, third and forth lines in the square parenthesis. Since the factor qip is obtained, in all 
other places of this formula qi can be set to zero. Then the (3/0) current contribution assumes the form 


3aH^o,t) = 
^ ^ ^ 
(3/0) 


dull 


duj2 


dui: 


4 

e 9s 

‘ 2S 


E 


i^qiLUi '^q2UJ2 V^qs^a 


p2(q3+qi+q2)-ro-'i(cJi+a;2+u;3)i 


E 

Zk 


qiq2q3 

(Zkllfia, e“*('l3+q2)-r|^|;^ k qg + q2) 


-9l/3 


Elk - Ei,k+q 3 +ci 2 + ^(wi -b a ;2 -b wa) -b ih'y Huji -b ih'y 


1 


dUlk - /bk-l-qa) 


ZD/k — ZD/^k+qa + h{uJi -b UJ 2 ) + ih'y 
1 


dk/3 

ZZ(/i,k+q3 ~ /i.k+qa+qa) 


ZD/,k-i-q 3 ~ ZD/^k+qa-t-qa + h{uji “b UJ 2 ) + ih'y 


dkp 


(B6) 


Now one sees that the terms with the Fermi-function differences in the square parenthesis are proportional to the 
wavevector q 2 times the second derivative of the Fermi functions. Omitting q 2 in the rest of the formula we get 

/ OO nOO nOO 

duJi / duJ2 / duJ3^ 

-00 J —00 J —00 

(3/0) 

(Zk|{{;a,e~^q=^ ‘~}+|Z,k-bqa) _ -qip _ -^27 d'^jfik - fiM+gs) 

ZD/k - ZD/.k+qa-b Zi/wi-b W2-b wa)-b *Zi7 fiwi-b ZZiy Zi(a;i-b a;2)-b iZiy dkpdk^ 

Finally, expanding the last Fermi-functions difference at small qa and introducing the Fourier components of the 
electric field, instead of those of the electric potential, = —iqcjiqui, we get Eq. (l54ll . 



E rh rh rh „dq3+qi-l-q2)'ro-i(‘^i-l-“2-l-i^3)t 
V'qi tdl Vci20J2 Y<43<^3 

qiq2q3 
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2. The contribution (2/1) 


The contribution (2/1) consists of the terms containing two Kronecker symbols in the matrix elements in 
Collecting the corresponding terms we get 


/ CO nCO nCO 

dhJl / duJ2 / duJ3—^ ^ 

qiq2q3 

^ (rk|K, k + gg + gi + gz) / (-1)^+^' ^ 

■“ -E'ck ~ -E'z,k+q3+qi+q2 + ^(^1 + UJ 2 + UJ 3 ) + iH'y I 2 


o*(q3+qi+q2)-ro-'i(cJi+u;2+cJ3)i 


ZZ'k 


E 


5u"'5i"n"qip'n^ + ^ii"'g27??k+qj<5i"i' + gsi^lk+qa+qi'^'"'*"'^'"'' 


fl'k — fl" 


■k+qi 


l"l" 


El'k — -E'i"',k+qi+q 2 + d{uJl + UJ 2 ) + 


- E 


Sw"dliiil>>q2jr]1 + gi/3??k+q2'^'"'' + 9357k+q2+qi'^'"''"'^'"' 


El'k - £^z",k+qi + ?ia;i + i/iy 
/Z",k+q 2 ~ /Z"',k+q 2 +qi 


Z"Z' 


- E 




I k/g- -I r-k/fc _ J/^^k+q2 ~ ^Z^^Mc+q2+qi _ 

E'/'k ^(^1 H“ ^2) “t“ E'/",k+q2 E'/^^^,k+q2+qi “h 

“I" ^ll'" //''',k+q3 .k+qs+qi 

E'^^/.k+qa ^Z,k+q3+qi+q2 h(^UJi “h UJ2'} H” E^Z^^.k+qs E^Z'" .k+qs+qi “h iH^ 


+ 


E 


'dv’ivqzs'nk + Sw"q2jVk+q3^i"i' + difiVk+q3+ci2^i'"i"^i''i' 


fr 


,k+q3+q2 


-fl. 


k+q3+qi+q2 


Z"Z" 


E^Z^^.k+qa E^Z.k+qa+qi+q2 ^(wi H“ ^ 2 ) “1“ 


^Z^^7k+qa+q2 -^Z,k+q3+qi+q2 "t“ “1“ ^^7 


(B8) 


Consider different contributions in Eq. (IB8E In the first line in the square parenthesis we have the Fermi-functions 
difference /z'k — /z".k+qi- The prefactor resulting from the matrix elements has two terms containing i5z"z' and one 
term which does not contain (5;//;/. The first two correspond to the intra-band contributions, while the last one - 
to the inter-band contribution. The same structure can be seen in the second, third and fourth lines of Eq. (IB8E 
Combining all intra- and all inter-band terms together and removing one of the Kronecker symbols in the intra-band 
sums we get 


( 2 / 1 ) 
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^ y (/'k|{D^, e--fe+qi+q2) r|^|^^ k + qa + qi + qz) / (- 1 )^+^' ' 
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E 


Z^ZZ^^^ Z/27l7k-fqi T Z73(517k-t-q2+qi 


fl'k fl'.k+ci3 
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El'k — £'z"',k-|-qi-(-q 2 + fl{uJl + IJJ 2 ) + ifCf El^]^ — E;/^k+qi + ^1 + 


'5zz>>g277k + g3^1?k+q2+q/^"^' _ 
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fl",k+ci 2 /z",k+q 2 +qi 


-E 

I" 

+ E 
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/z",k-|-q 3 /z",k-|-q 3 -t-qi 
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/z'k fl"M-\-cii 


E^Z'k — E^Z'",k+qi+q2 + h{uJi + LO2) + E^Z^k — -^Z",k+qi + + ^^7 


- E 

l"l"' 
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Sii">qii3r]k+(i2^i"i' 


/z",k+q 2 /z"',k-|-q 2 -|-qi 


liijfi El'k — £'Z"',k+q2-|-qi + d(ujl -k LO2) + ih'y El;//^k+q2 ~ T'Z"',k+q2-(-qi + + *^7 

_ Z^ZZ'" qip 7 k-rq 3 _ /z",k+q3 /z"',k+q3-Sqi _ 

liqfii '®Z",k+q 3 ~ '®Z,k+q 3 +qi-|-q 2 + h{uil -k UJ2) + *^7 ^Z",k+q 3 ~ K^Z"',k+q 3 +qi + EOl -k i/17 
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E 




/j'",k+q 3 +q 2 /i,k+q 3 +qi+q 2 


-E';",k+q3 ^'i,k+q3+qi+q2 + h[uJl + UJ 2 ) + ifvj i?;'" ^k+q3+q2 -E'/,k+q3+qi+q2 + + *^7 J 


(B 9 ) 


Now one sees that each of the intra-band terms contains the Fermi-functions difference of the type (/(".k-t-qs ~ 
//",k+q3-i-qi) ~ —<lipdfi'i;if^+^^/dkp and all the inter-band terms have the prefactor qip originating from the matrix 
elements. One can therefore set qi = 0 in the rest of the formula, which simplifies the second energy-denominators, 
e.g., ifi",k4-q2 ~ ^^i",k+q2+qi + ^1 + ifi'ck+qa ” ^^i",k+q2 + + *^7 = ^1 + Then, taking the sums over 

I" and I"' in the inter-band contributions we present the result in the form 


^ ^ > 

( 2 / 1 ) 


/ OO pOO nOO £>^ n 

duji / duj2 / ^ 
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Sii"q2 ^vZ+ q3svl+q2^i"i' 


-gl/3 

fvjJi -b itvy 


E 


9(/z'k - /i'bk+qs) 


F'Z'k — £'Z",k-i-q 2 + ^(wi + 1112) + id'y 
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^W'QSS'dk. 5 (/z".k+q 3 //,k+q 3 -|-q 2 ) 


■Jp i/;",k-|-q3 - ■E/,k-(-q3-fq2 + ^(^1 + ^ 2 ) + iflj dkp 

^ (/C^,(k)-/C'!,(k + q 2 )) 


J =A 


gl /3 


El’ln — ifz,k+q 2 + + W2) + *^7 

1 

T'Zyk+qs ~ £'Z,k+q3-|-q2 + ^(wi -b W 2 ) + *^7 
where the function /Cpj(k) is defined as 


(Ez(^ + ‘la) ~ Ez(^ + q 2 + qs)) 


J =B- 


Ez(k)=i 7 ^ 


/z'k — /zi< 


(BIO) 


(Bll) 


■ ifz/k - ifzk + + *^7 

The terms in the first square parenthesis in (IBlOll denoted as A are transformed as follows. First, we take the sums 
over I" and set q2 = 0 or qa = 0 in the terms which have already the prefactors <727 or q^s'- 


A = 


927 gk 


d(fi'k - /zk) 


927 gk+ 


qs 


^(/Z',k+q 3 /z,k+q 3 ) 


F'z'k - £^zk + ^(wi-b a; 2 )-b 1/17 ^k/s Ei>^k+qa - Ei^k+q^ + + UJ 2 ) + ifij dkp 

93517k 5(/zk ~ /z,k-|-q 2 ) 935 l?k+q 2 5(/z'k ~ /z',k 4 -q 2 ) 


Elk — El^k+q 2 + + OJ 2 ) + *^7 


dk^ 


El'k — £'Zbk+q2 + d(uJi -b UJ2) + ikrf 


dk/3 


(B 12 ) 


Now consider the first line in Eq. (|B 12 I) which has already the prefactor 927- One sees that the second term there 
is the same function as the first term but with the argument k -b qa. It can be expanded at qa —> 0 which gives the 
prefactor q^s- In the second line the differences of the Fermi functions can be expanded at q2 —>■ 0 ; in the rest of the 
second line we can then set q2 = 0 . So we get: 


Asi - 


927935- 


1 

9 k 


i9(/z'k — /zk) 


A 

dks \Ei'k — Elk + H^Ji -b C02) + ih'j dkp 


927935 


9 ^ 


d^iflk — fl'k) 


h{uii -b a;2 -b *7) dk^dk^ 


(B 13 ) 


The terms in the second square parenthesis in (IBIOI) denoted as B are transformed as follows. First, we expand the 
difference of the /C-functions at small qa and set qa = 0 in the rest of the formula: 


B 


I 


-927 


\ T^z'k — Elk + ^(wi + 0J2) + ih'j dkj 

_1_ d 

El', k+qa - El, k+qa + + ^2) + 1^7 9k. 


^ (Ez(k)) 

;(Ez(k^ 


q 3 


(BI 4 ) 
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Now we see that the expression in the second line is the same function as in the first line but with the argument 
shifted by qs- Expanding (IB 14 I) at qa —>■ 0 we then get 

„ _ _^( _ \ _^ f p fi'k - fik 

927<735 ^ + a;2) + ih'y dk^ - Ei]^ + fiwi + ihj 

One sees that both the A and the B terms differ from zero only if V = I, i.e. A,B (x 5 i,i. Substituting the terms A 
and B in (IBIOI) . taking the sum over V and introducing the Fourier components of the electric field as was described 
in Section ED we get the expression (ESI). 



3. The contribution (1/2) 


The contribution ( 1 / 2 ) consists of the terms containing one Kronecker symbol in the matrix elements in 
corresponding terms are 


/ OO nOO nOO Q 

dujl I dui2 I duJ3 'y ^ 

-OO J —OO J —OO „ „ „ 


pi(q 3 +qi+q 2 )'ro-i((^l+ 1 ^ 2 +1.^3)* 


( 1 / 2 ) 


y (l'k|{D^, e--(q3+qi+q2)-r|^|^^ k + qa + qi + qz) (-1)^+^' 

Ein^ — Ei^k-fq3-fqi+q2 + ^(<^1 + W2 + wa) + ih'y 4 


ZZ'k 


Sii'''qipq 2 jVk'nl 


/z'k — //"k 


^'/'k ~ -E'/"'k + h{u!i + IJJ2) + ih'y E;'k ~ + hioi + ih'y 


E 

i"i"' 

y - 

-^''k - ■Ez"yk+q 2 + + 1 ^ 2 ) + ih'y El'k - Ei'^w + huJi+ ih'y 

A 


gi/ 3935 l?k-rq 2 '^'"''"'^k 


/z'k — /z"k 


g 27 g 3 ( 5 l 7 k+qi ^k+qi _ /z'k — /z",k-rqi _ 

-^i'k - -Ez-'-.k+qi + h{LOi + W2) + ih'y Ej/k - -Ei//,k+qi +huJi+ ih'y 


E 


- E 


Su'"qipq 2 jVk'nk 


/z"k — fl" 


-^''k - El'i'k + h{uji + UJ2) + ihy Ei»i^ - Ei>»u + hwi + ihy 


- E 


g27g3^1?k-fqi hk _ /Z''k — //''',k-Tqi _ 

liTjfri -^^'k ~ -E'P'yk+qi + h{uji + LU2) + ihy Eii'k — i?z"',k-i-qi + huji + ihy 

qipqAhi+^^hi+^Jv'i' 


- E 


/z",k-|-q2 - /z"'.k+q2 


-^i'k - Ei"'.yi+^2 + +1^2) + ihy k+q2 - -E////,k+q2 + huji + ihy 

fl"W — fl"'k 


SwqipqssVkVi 


- y _ 

iTijf,, Ei"w - -E/,k+q2 + ^(wi + W2) + ihy Ei^i^ - E;///k + hwi + ihy 

_ g27g3^1?k4.q^</z"'Z"l?k _ /z''k ~ /z"',k-rqi _ 

Evk — E/^k+qi + h{uJi + UI2) + ihy Ei"\^ — E;'"^k-fqi + Eoi + ihy 


- E 


gl/ 3 g 27 » 7 k+q 3 l 7 k+q 3 '^i"*' 


/Z",k+q 3 /Z"'.k+q 3 


-E'Z",k+q3 ~ -E'z,k+q3 + h{uji + 0^2) + ihy i?i"^k+q3 ~ -£'Z"yk+q3 + EjJi + ihy 

_<^ZZ"'g 27 g 3 i 5 l?kl?k_~ /z.k+qi_ 

Ei"w - Ez,k+qi + h{uJi + LO2) + ihy Ej/'-k - -Ez,k-i-qi +hhji+ ihy 


gl^g 3 ( 5 ??k-|-q 2 ^i"'^"llk _ /z'",k+q2 ~ /z,k+q2 _ 

Evk - -Ez,k+q2 + ^(wi + W2) + ihy £;;///^k+q2 - -E^z.k+q2 + fiwi + ihy 


E 


_ ^27 ^k+qa ^k+qa _ /z^^^k+qa //,k+q 3 _ 

£;///,k+qg - ^Z,k+q3 + ^(^1 + ^ 2 ) + £;////,k+qa “ ^Z,k+q3 + ^Wl + ^^7 


E 


. The 


(B 16 ) 
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where we have set qi = 0, q 2 = 0 or qa = 0 in the rest of the terms if they already contained the corresponding 
q-prefactors. In each of the twelve terms in the square parenthesis in (IB16I) we have to sum over V" and I". This can 
be done as follows. Consider, for example, the first line in the square parenthesis in (IB16|) . The sum over V" is taken 
due to the Kronecker symbol; I"' is then replaced by L To take the sum over I" we notice that I" takes two values, 
for example, /' and I'. As seen from the Fermi-functions difference, /i'k — /;"k, the first line vanishes if I" = I'. Hence, 
I" must be equal V. Similarly, one can simplify eight terms out of twelve. 

Now consider the third line. Here I" = I' due to the Kronecker delta. The Fermi-functions difference can then be 
expanded at qi —>■ 0 and qi can be set to zero in the rest of the formula. Using this way one can simplify the other 
four terms. Then we obtain: 


/ CO pCO p 

dwi / du )2 / 

-CO j —CO J — 


( 1 / 2 ) 


oo 4 


i:' 

qiq2q3 


H^c[ 2^2 V^qa^^a '• 


^i(q3+qi+q2)-ro-2(a;i+a;2+cc;3)i 


E 


(rk|{D„, e-*('i3+qi+q2)-r}^|;^ k -h qg -f qi + q2) (-1)'+'' 


El'\i — Ai^k+q 3 +qi-|-q 2 + ^(^1 + W 2 “k UJ 3 ) “k itl'y 4 


gi/3'?27 


(3 7 


fl'k - fv 


El'k - -k h{uii + U 2 ) + ih'^ Aj/k - Upk -k hui + ih'j 




fl'.k+qs fl',k-\-q 3 


El'^k+qa — Aj^k+qa + h{uJl + UJ 2 ) + iEy Aj'^k+qa ~ ^i',k+q 3 + ^1 + 

flk - flk 


'?1/3<?27 I 


p 7 

€Vk 


A;'k — Ajk + h(uj\ -k W 2 ) + ifi'y Pjj — A/k + Eul “k ih'y 


^k-l-qa^k-l-qa 


fLk-'rqa flM+qs 


El'M+qa ~ ^Ck+qa + + UJ2) + iEj p — Ai^k-l-qa + EjJi + iH'y 


+ 


qimsv^+c 


Vk- 


fl'k — fl'] 


El'k - Ap,k+q2 + ^(^1 + ^^ 2 ) + ih'Yy Ei'i^ - Ef,^ + hui+ ihj 

/i',k+q 2 -fv ,k+q 2 \ 


^k+ 


q2 


Az/^k+q2 - ■E^/'.k+q2 + Swi -k iH'y 


qimsdk 


Ei^ — Ui^k-i-qa + + W 2 ) + ^ik ~ Eii^ + huji + iH'y 

fl,k+qa /kk-|-q 2 \ 


flk — flk 


7k+ 


^r,k+q 2 - El,k+qa + Swi + ih-y 


-qifiq2jq3svivl 


HuJi + ih'f 


E 


difvk — fv'k) 


El'k - Ei"k + fi{ui + UJ 2 ) + ih'f dkp 


-E 


d{fi"k — flk 


-jjf Ei"k - Elk + fi{u}i -k a;2) -k ih-y dkp 


(B17) 


in (|B17I) we have also rearranged terms and collected them in groups according to their g-prefactors. 

Now one sees that the function in the second line differs from the function in the first line by only the shifted (by 
qa) argument (k —k -k qa). It can therefore be expanded in qa and simplified. The same can be done in all other 
lines except the two last ones which still contain the sums over These sums can now be also easily taken since in 
the last line I" should be equal to /, while in the last but one line - to V. Making all these transformations we get 


( 1 / 2 ) 


dull 


duj 2 


dui: 


4 

e 9s 
' 2S 


E < 

qiq2q3 


«^qia;i V^q2CL!2 V^q3(^3 


o*(q3+qi+q2)-ro—'i(a;i+a;2+a;3)i 


(/'kllDc,, e-^(q^+q^+q^)'P+|/, k -k qa + qi + q 2 ) (-1)^+^' 
El'k — Ai^k+q3+qi-|-q2 + ^(^1 "k W 2 "k Wa) + ifl'J 4 


ll'ic 


X 
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„ „ „ ^ 

^ hihl 

//'k - fl'k \ 

91/3927935^^^ 

\^Ei’k — Elk + h{uji + W 2 ) -b ih'j Eiik 

— Epi^ + hijji + ih'j j 


91/3927935 


d ( _?49k_ fik - fik 

dks I i^i'k — Elk + ^(wi + W 2 ) + ih'j Ej]^ — Elk + Eoi + ih'y 


- 91/3927 935 %+q 2 


d 


Ei'k - Ep k+q^ + ft(wi + 0 / 2 ) + ih'j dk^ [^^Ei^k - Epy. + Hui + ih'j 


/i'k - fv 


- 91 / 3927935 % 


d 


Elk ~ ^'/.k+qa -I- h{uJi + W 2 ) -I- ih'j dk-y ^^ Ejy. — Elk + hui + ihj^ 

-9l/3927935l?k??k [ _ I _ ^(//'k - fl'k) 

huJi+ih'j \Ellk — Ef,^^ + h{uJl+UJ 2 )+ih'J dkp 

_1_ djfik-fik) ' 

Elk — Elk E h[uii + W 2 ) -l- ih'j dkp 


fik - /ik 


(B18) 


Now we have the required three g-prefactors 91^927935 and can set qi = q 2 = qa = 0 in the rest of the formula. Then 
it can be rewritten as 


j^'>{ro,t) = duji duj 2 I dw. 


00 • 4 

le^gs 


( 1 / 2 ) 


-00 «/ —00 


2S 


V E^ E' 

/ V qi(^i q2CL!2 < 


'<5 p^(q3+qi+q2)-ro-'i(a;i+CLJ2+u;3)i 

C130J3^ 


qiq2q3 




{l'MVa\lk) 


(- 1 ) 


l+l' 


^ El'k - Elk + h(uJi + W 2 -I- W 3 ) + ih'j 2 


\ ^ 

^ hihl 

1 

flk - flk \ 

dks 

y^Ei'k — Elk “b h[uji + W 2 ) -b ih'j 

^Ei^k — Ej,^ + hijJi + ih'j 

Efk — Elk + hui + ih'j J 


-b 


€ 


d 


T)^ - 

£^i'k — Epi^ + h{(jji + 0/2) + ih'j dky I Ei^k — £^pk + ^1 + ihl 


fi'k - h 


hi 


d 


Elk ~ Ei,k + h{uJi + UJ2) + ih'j dkj [^^ Eiy. — Eik + fiwi + ih'j^ 

1 dif,k - /Pk) 


flk ~ //k 


9^/k 


hull -b ih'j \ Ei^k — Epk -b h{u!i + UJ 2 ) + ih'j dkp 

_1_ djfik-fik) ' 

Elk ~ Elk + h{uji + 0 / 2 ) -b ih'j dkp 


(B19) 


Now we can take the sum over V. One sees that, if /' = / all the expression in the square parenthesis vanishes. 
Hence I' should be equal to I and we get 


;l3)(ro,t) = / dui, I duj2 j_ da33^<(ro)£;3.(roX(ro)e-*(-^+-=+-=')‘ 


( 1 / 2 ) 


' —00 —00 


E 

/k 

(/k|i)a|Zk) 

\ ^ f 

hihi 


flk - flk \ 

h{uJi -b 0/2 -b 0 / 3 ) + ih'j 

dks y 

h{uji + UJ 2 ) 

+ ih'j Elk 

— Eji^ + fujJi + ih'j j 


hi 

d 


flk — flk 

\ 


-b 


Elk — Ejy. + h{uji + UJ2) + ih'j dky I Elk — Ejy. + hioi + ih'j ^ 

1 difik-fik)' 


hivi 


hui + ih'j \ Elk — £^rk + hiuji + W 2 ) -b ih'j dkp 


(B20) 


Since the electric field, and hence the current, does not depend on the coordinate, the expression (IB20I) gives the 
formula ((551) . 
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4. The contribution (0/3) 

The (0/3) contribution is the only one term in Eq. (IB3|) which contains no Kronecker symbols. This does not 
allow one to immediately take some of the sums. On the other hand, the (0/3) contribution has all three g-prefactors, 
Qil 3 <l 2 -'fQ 3 Sj so that one can set qi = q 2 = qs = 0 in the rest of the formula and to replace the Fourier components of 
the potential by those of the electric field. So we get straightforwardly: 

/ OO j-OO nOO . 4 

-cx) J —00 J —00 


( 0 / 3 ) 


E' 




Ei'k - Elk + + a;2 + wa) + ih'y 


E 

i"i"' 


1 




fl'k — fl"'V 


El'k - Ei»k + + W2) + Ei>k - Ei>»k + + ifij 


- 


1 


fl"k — fl'"l 


l"V‘ 


El’k - Ei^k + hiioi + W 2 ) + ih'j El,Ik - Ei,„k + Eoi + ih'j 


- E 


1 


//"k — fl"'k 


l"l'‘ 


E 


Ei''k ~ Elk + + W2) + iEy Ei,,k — Eimk + huji + ih'j 


1 


fl"'k — flk 


Ei"k - Elk + + W 2 ) + ih-f Ei,„k - Elk + fiwi + ifij 


(B21) 


Now we take a sum over I'". In the first line in the squared parenthesis I'" should be equal I', in the second line - Z", 
etc. We get 


^ = 
( 0 / 3 ) 


dcoi I duj 2 I 


— 00 J —00 




(Z'k|{)„|Zk) 




Eiik — Elk + Zi(wi + W2 + oJs) + ih'y 

1 fl'k — fl'k 


E 




£^//k — + h{uji + CJ 2 ) + + ^^7 


- E 


fi"k - fi'i 


- E 


■jjf Ei,k - Ei„k + Zi(a;i + 0 ^ 2 ) + ih'j Ei„k - Ei„k + fiwi + iZiy 

_ \ _ /i"k - fl"k _ 

7 Ev'k - Elk + h{uji + W 2 ) + ih'y Eiiik - Ei„k + Eoi + ih-f 

1 flk - flk 


E 


Ei"k — Elk + + <^ 2 ) + ifi'l Elk ~ Elk + Eoi + ihj 


(B22) 


Now we replace, in the first line in the squared parenthesis, the summation index I" by I"; this does not change the 
result. Then the (0/3) current is written as 

o4. 


/ OO POO poo ■ 4 

dwi / d0J2 / doj.^E^^EZ^El 

-OO j —OO j —OO 

(/'klucIZk) 


1(5 —i(l^i+lil 2 +<.iJ 3 )t 

bjp 


( 0 / 3 ) 


p M -L -mr^~^E'vtvkVk 

Ei'u_ — Eiiti + h{uJi UJ2 0 J 3 ) + ^^7 


E- 


fl'k — fl'] 


fl"k — fl" 


El'k — Ei„k + fi{uJi + UJ2) + ifij \Ei,k — Epk + fnoi + ih'j Ei,,k — Ei„k + huji + ih'y 


E 


1 


fi"k - fl" 


flk ~ flh 


-jjf Ei"k — Elk + Zi(wi + a;2) + ih'y \Ei"k — Ei„k + fpJi + ih'y Ejk — Eik + Ziwi + ihj 


(B23) 


X 




























31 


Now one sees that I" should be equal to V in the first line and I in the second line; otherwise, the corresponding terms 
vanish. This gives 


/ OO pOO pOO • 4 

/ cL;, / 

-OO j —OO j —OQ 

(0/3) 


X 


X 


:E' 


{l'k\vc.\lk) 


h{uji + UJ 2 ) + ih'y ^ Ei>k - + h{uji +U2+ los) + 




fi'k - fi>: 


fi'k - fi'k 


Ei'-k - Ej,y. + HuJi + ih'y Ej,^ - Ei>\^ + huji + ih'j 
/ik — /tk /tk ~ fik 


Elk — Efy. + hwi + 1/17 Ej]^ — Elk + Eoi + ih'j 


(B24) 


Finally, in this formula, if F = / the term in the squared parenthesis vanishes. Thus V should be equal 1. Substituting 
V = I we get, after some transformations, Eq. (1571). 


Appendix C: Calculation of conductivity 

In this Section we calculate the contributions to the third-order conductivity defined by Eqs. (l58l) and (l54ll - (1571) . 


1. The contribution (3/0) 


Following the definitions (l58l) and (|54)) we present the (3/0) contribution to the third-order conductivity in the form 


^a/375(wi>‘^2,W3) = 


w, 


(3/0) 

(y/3'yS 


(ill + ^ 1^2 “t“ + ?r)(rii + ^2 xr)(r2i + ^r) 


(Cl) 


(3/0) 


where the factor is 


W, 


(3/0) _ 


aB-yS IT’S c 

Ik 


dkpdkjdks 


(C2) 


Substituting in (IC2I) the matrix element of the velocity (1251) and integrating by part we rewrite the function 
in the form 


>,y(3/o) _ is'^ds df{Eik) dEik Elk 

~ hE%S ^ dE dks dkc^dkpdk^' 


(C3) 


Due to the factor df{Eik)/dE only the vicinity of the Fermi level contributes to the integral. Therefore the integration 
over the whole Brillouin zone in (IC3I) is reduced to the integration over the vicinity of two Dirac points. The both 
Dirac points contribute equally to the integral which leads to an additional, valley degeneracy factor = 2. Near 
the Dirac points we have Eik = {—lyhvpk, Eq. (IT^ . The differentiation gives 


dEik 

dks 


{-lyhvF^, 


(C4) 


d^Eik 

-p- + (“) 

Substituting (IC4|) and (IC5I) back to Eq. (IC3I) and calculating the integral leads to the following result (at T <C |/r|): 

Here the angular integration in the k-space [k = fc(cos </>, sin (/>)] has been performed using Eqs. (IA3I) and (IA4I) . With 
the help of the definitions (ISTl) , (1551) and (IMl) we get Eq. (I55|) . 
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2. The contribution (2/1) 


According to the definition (l 55 l) the (2/l)-conductivity has three contributions, 


( 2 / 1 ) 


(2/1),a 

■E- 


(2/l).fc 

(rk|z;„|;k) 


(2/l).c 


• 4 

2S ^ — iJjk + h{ijJi + UJ2+ wa) + ih'y 1 + uj 2 + ij) 


s d^ifik - /ik) 


Vk 


dkpdk^ 


+ 


d 


O' 


difik - flu) 


h{uji + i"/) dks I — Eiu + h{uii + UJ 2 ) + ifi'j dkp 


d 

dks 


1 


d 


Ely. — Elu + h{uji EUJ2) + *^7 dk.y Ely. — Eiu + ftwi + ih'y 


fik — fiu 


(C7) 


corresponding to the three lines in Eq. (IC7l) . Consider them one after another. 
The first-line term in Eq. (EH) can be written as 




1 


• 4 

le gs 


E 


(Zk|i)a|Zk) 


sd'^ifik- fik) 
Vk- 


4:OiOi2 2EpS it'jk — Eiu -f h{u}\ E UJ 2 + W 3 + * 7 ) dkpdk. 


(C 8 ) 


(2/1),a 

Substituting the inter-band matrix element of the velocity (|26|) and taking the sum over I we obtain 

1 ie'^Qs ( 2.Eu 2,Eu 


4175‘^ 2 , 013 ) = 


4O1O12 AhElS 


E 


2 Ek -I- ?i(a;i -|- a;2 + W3 + *7) 2 £’k — ?i(wi -I- u;2 + W3 + *7) 


(2/1),a 




(C9) 


dkpdk.y 

where Ek = E 2 u- The Eermi-functions difference in (IC9I) . (/ik — /2k), is proportional to Q{k—kp) at T = 0, therefore 
the integrand in (IC9I) contains the derivative of the delta-function 5'{k — kp), only the vicinity of the two Dirac points 
contributes to the integral, and the integration over dk can be easily done. Performing the angular integration with 
the help of formulas (jASp and (|A 4 p we get 


414 ^ 2 ,^ 3 ) = "^4 


4C?iC?i2 


( 2 / 1 ),' 


^aS ^13') 


1 


1 


“t“Oi23 ^^a5^/37 


1 — O 123 1 + O 123 

1 


This gives the first (a) contribution to D 2 ,D 3 ). 

The second-line term in Eq. (IC7I) has the form 

414(‘^T‘^2,W3) = 


ie^'^gs 


:E' 


(1+0123)^ (1 — 0123)^ 


{Eju - Eiu)g^ 


(CIO) 


(2/l).fc 


SOlhEpS ^ Ely. - Eiu + ft(wi -I- a ;2 + wa -I- 17 ) 


7 

Vu 


d{fik - fik) 


Calculating it similarly we get 

445 (^ 1 : ^^ 2 ,^ 3 ) = 4 ^^ 


(2/l).b 


A _ 

dks \ Efi^ — Eiu + h{uji -I- a ;2 + * 7 ) dkp 




(Cll) 


1 


“I" I ^a.^^'yS 


^07 


(1 + Oi2)(l + O123) 

1 1 


(1 — O123) / ~ Ol2)(l — O123) j 


(C12) 
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This gives the second (6) contribution to ri 2 , 

The third-line term in Eq. (IC7l) . 


-(3) / N ie'^gs 


(^Fk - Eiv,)v^ 


( 2 / 1 ).' 


dfiS* — Ei\^ -\- h{uji UJ2 “1“ W3) -|- ifvy 


dks 


d 


Ef^ — Eii^ + h{uji -I- W2) -I- ih'^ dkj Efy. — Ei]^ + hcui + ih'y 


/ik ~ /ik 


can be reduced, after a long integration, to the form 


~ 3) , \ (3) '■ 

</375(w1’‘^2,W3) = CTq 7 
'»_ _✓ ^ 


( 2 / 1 ).' 


{SapSjS — Aa,fl 75 / 4 ) ^Jl(Ol, Ol2, O123) — Jl( —Ol, — O12, —Ol 23 )^ 


{Sa’ySpS — ^aP’ys/'^) (yJ2{Ol,Oi2,Oi2'i) — j2{ — Oi, —O12, —Oi 23 )^ 


(C13) 


(CM) 


where the integrals Ji{a, b, c) and 1/2 (a, b, c) are defined in ( 1551 ) - ( 170 )) . The formula (IC14I) gives the last (c) contribution 
to O2) ^^ 3 )- Combining all three contributions we get after some transformations Eq. ( 1001 ) . 


3. The contribution (1/2) 

According to the definition (ISOl) the (l/2)-conductivity also has three contributions. 


^a/ 375 (‘^b ^^2, W3) = a-),^^.^5(u;i, W2, wa) -f W2, wa) -f W2, wa) 


A3) 


A3) 


A3) 


( 1 / 2 ) 


(l/2).o 

1 ^e gs 
O123 4 :EpS 


(l/2).& 


^(^k|fi„|Fk) 


Zk 


(l/2),c 
d f f3 J 


/zk — /zk 


h{uji + UJ2) + ih'y dks I Eiy, - Ej^ + hwi + ih'y 


hi 


d 


Elk — Ej y. + h{uJi -I- a; 2 ) -f ih'y dk^ Eik — -|- fiwi -|- ih'y 

1 difik-fik)' 


fik - fik 


hZhi 


huji -I- ih'y Elk — Ei^ + h{uJi -I- W2) -I- ih'y dkp 
corresponding to the three lines in Eq. (IC15E The first one (a) has the form 
-(3) , s 1 ie-'^gs sr^ {lk\va\lk) d ( p ^ 


(l/2).a 


J, O123 4£’i7’5' ^ h{uji -I- UJ2) + ih'y dks 


hkVk 


fik — fik 


Elk - Ely. -I- hioi + ih'y 


Calculating it like in the previous cases we get: 


^aj75(wi,W2,W3) = (7^^ 


8O123O12O1 \ 2O1 1 — Oi 


' lni±^-l . 


(l/2).a 


The second (6) contribution is 

1 


~(3) , x_ 1 ie.-gs ^ {lk\va\lk)g^ d 

'^0!/375('^1’^2j W 3) n _ T?.. T?- I k!,,. i i Ak^, au I hk 


fik — fik 


✓ O123 4 :EpS ^ Elk — El y. -I- h{uJi -I- 0:2) -I- ih'y dk.y I ^ Eik — Eyy + hui + ih'y J 


(C15) 


(C16) 


(C17) 


(C18) 


(l/2).b 
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Calculating it we get 


W 3 ) = 40^23 I ~ ^“/57i 5/4) (^J3(Ol, O 12 ) + J3(-Ol, -O 12 )) 

(1/2),b 

— ^SajS/SS — Aa0js/'i'j (^J4(Oi,Oi2) + J4( — Oi , — Oi2)^ 

where the integrals Jsia, b) and J 4 (a, b) are defined in Eqs. (1711) - (1771) . 

The last (c) contribution corresponding to the third line in (IC15I) can be written as 

~(3) , 1 ie'^gs {lk\vc\lk)r]lrjl '9(/ik -/fk) 

OlOl23 4£;2 5^£;^k-£;^k + ^(wl+W2) + ^?i7 dkp 

(1/2),C 

After the integration we obtain 

^i%si^l^^2,i^3\ = 40^2301 (tT^ ■ " ^«/37^/4) ■ 

(l/2).c 


Combining all three contributions we get after some transformations Eq. ([67|. 
The last, (0/3), contribution is calculated similarly. 


(C19) 


(C20) 


(C21) 


Appendix D: Special cases for the integrals Jn 

The formula (17^ - (ITSl) provide results for the integrals Ji{a, b, c) - b) at all values of the arguments a, b and 
c but in the “degenerate” cases a = b, a = c, etc. it is useful to have their simplified explicit expressions. Here they 
are. 


1. Integral (a, &, c) 


In the special “degenerate” cases we have: 


Ji(a, a,c) = -^ + 


c — 2a 


4c — a 


1 


1 


a^c a2(a — c)^(a + 1) c{a — c)^{c+l) 2(a — c)^(a + 1)^ (a —c)2(c+l)^ 

—3a^ + 12a^c — 8ac^ + 2c^ , , , 2a — 5c , , ., 

-ln(a + 1) - ^ ln(c + 1), 


a^(a — c)^ 


c(a — cY 


(Dl) 


Ji(a, 6, a) = -^ + 


2a — b 


^b a(a —6)^(a + l) (a —6)^(&+l) 2a(a — 6)^(a + 1)^ 3(a —6)(a+l)^ 


—a^ + 2a^6 — 3a&^ + 6^ — Sab + Sb^ 

- ln(a + 1)- \n{b + 1), 


a3(a-6)4 


62(a — 5)^ 


(D2) 


Ji{a,b,b) = -75 + 


a — 46 


46-a 


ab"^ 6(a —6)3(6+1) 26(a — 6)^(6 + 1)^ (a —6)(6+l)3 

—a^ — Sab + 6^ a^ — Sa^b + 2ab^ + 36^ 

■ ln(a + 1)-pr---ln(6 + 1), 


a^{a — 6)"^ 


63(a-6)4 


(D3) 


^ ^ 1 1 
Ji(a, a,a) = — + 


3(a + l) 4(a+l)^ 3a(a +1)3 


1 2 , , 

- ^Ma + !)■ 


(D4) 
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2. Integral ^2(ffl, 6, c) 

For the second integral in the special cases we have: 


J 2 {a,a,c) = 


2a — c 


2a(a — c)(a + 1)2 0^(0 — c)2(a + 1) a^{a — c)^ 


3a^ — 3ac + 1 

ln(a + 1) + ^ ln(c + 1), 


c(a — c)^ 


(D5) 


J 2 {a, b, a) = -- 


,(a —6)2(a + l) 6(a —5)2(6+1) a‘^{a — b)^ 


b 3qj ,, , 35 a ^ \ 

ln(a + l)-T3-^-;^ln(5+l), (D6) 


52(0 — 5)3 


^ 2b—a 1 1 1 / + 3a5 — 352 

^2(0, 5, 5) = - - —- - - -ln(a + 1)- ^ -ln(5 + 1), (D7) 


52(0 — 5)2(5 + 1) 25(a — 5)(5 + 1)2 a{a — b)^ 


b^{a — 5)3 


J2{a,a,a) = -- 


1 1 , / 

+ ^ln(a + l). 


i3(a + l) 2a2(a + l)2 3a(a +1)3 a'^ 


(D8) 


3. Integrals J 3 (a, 5) — J 5 (a, 5) 


For these integrals the number of special cases is small and we get: 

2 11 
7/3(0, o) =-7 ln(l + o) H—^ + 


774 ( 0 , 0 ) = - 


o2 o2(l+o)’ 


H—7 ln(l + o) 


2a(l + o)2 o2(l + o) o3 


(D9) 

(DIO) 


4 2 

175(0,0) = -775(0,-0) = -7 + 


Ai (^-°) 
o2 ' o2(l — o2) o3 ( 1 + 0 ) 


(Dll) 
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